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Preface 


In  late  1976,  a  study  to  produce  a  wave  climate  for  U.  S.  coastal 
waters  was  initiated  at  the  U.  S.  Army  Engineer  Waterways  Experiment 
Station  (WES) .  The  Wave  Information  Study  (WIS)  was  authorized  by  the 
Office,  Chief  of  Engineers,  U.  S.  Army,  as  a  part  of  the  Field  Data 
Collection  Program  which  is  managed  by  the  U.  S.  Army  Coastal  Engineer¬ 
ing  Research  Center.  The  U.  S.  Army  Engineer  Division,  South  Atlantic, 
and  the  U.  S.  Army  Engineer  Division,  New  England,  also  authorized 
funds  during  the  initial  year  of  this  study  (FY  1978)  to  expedite  exe¬ 
cution  of  the  Atlantic  coast  portion  of  this  program. 

This  report,  the  twelfth  in  a  series,  presents  a  simplified 
version  of  the  wave  hindcast  model  developed  to  calculate  the  wave  in¬ 
formation.  The  study  was  conducted  in  the  Hydraulics  Laboratory  under 
the  direction  of  Mr.  H.  B.  Simmons,  Chief  of  the  Hydraulics  Laboratory, 
Dr.  R.  W.  Whalin,  Chief  of  the  Wave  Dynamics  Division,  and  Mr.  C.  E. 
Chatham,  Jr. ,  Chief  of  the  Wave  Processes  Branch.  This  report  was 
prepared  by  Dr.  D.  T.  Resio  and  Mrs.  B.  A.  Tracy  with  assistance  from 
Mr.  W.  D.  Corson,  Mrs.  D.  S.  Ragsdale,  and  Mr.  R.  E.  Jensen.  Dr.  C.  L. 
Vincent  provided  valuable  suggestions. 

Commanders  and  Directors  of  WES  during  the  conduct  of  the  study 
and  the  preparation  and  publication  of  this  report  were  COL  John  L. 
Cannon,  CE,  COL  Nelson  P.  Conover,  CE,  and  COL  Tilford  C.  Creel,  CE. 
Technical  Director  was  Mr.  F.  R.  Brown. 


Summary 


In  1976,  the  U.  S.  Army  Engineer  Waterways  Experiment  Station 
began  a  study  to  produce  a  wave  climate  for  U.  S.  coastal  waters.  This 
climatological  information  is  produced  by  numerical  simulation  of  wave 
growth,  propagation,  and  decay  driven  by  historical  wind  fields.  It  is 
imperative,  if  such  an  approach  is  to  be  used  for  applications  of  sig¬ 
nificant  economic  consequences,  that  the  entire  set  of  input  data,  all 
numerical  techniques,  and  all  general  assumptions  be  thoroughly  investi¬ 
gated  and  documented  to  determine  the  types  and  magnitudes  of  errors 
intrinsic  to  their  use. 

There  are  four  basic  steps  in  the  calculation  of  waves  from  past 
meteorological  data.  First,  pressure  data  must  be  assimilated  into  a 
pressure  field  that  depicts  all  important  synoptic  weather  features. 
Gradients  of  pressure  in  time  and  space,  along  with  certain  thermal 
characteristics  of  the  planetary  boundary  layer,  are  then  used  to  con¬ 
struct  an  estimate  of  a  quasi-geostrophic  wind  speed  and  direction  at 
some  level  where  it  is  assumed  that  the  frictional  effects  of  the  ocean 
surface  on  the  atmosphere  are  negligible.  Next,  an  analysis  of  the  ver¬ 
tical  variation  of  the  wind  in  the  planetary  boundary  layer  is  used  to 
reduce  this  wind  to  a  common  19.5-m  level.  Finally,  these  surface  winds 
are  input  into  a  numerical  wave  model  to  simulate  wave  generation,  prop¬ 
agation,  and  decay. 

If  any  one  of  the  above  steps  contributes  significant  bias  (on  a 
geographical  basis,  seasonally  or  overall),  it  can  introduce  errors  into 
the  results  that  are  difficult  or  impossible  to  remove.  Similarly,  if 
any  step  contains  a  large  random  error,  certain  statistics  (such  as 
duration  curves,  extremes,  and  conditional  probabilities)  can  be  seri¬ 
ously  affected.  Thus,  each  step  must  be  checked  independently  where 
possible.  This  serves  to  substantiate  the  merit  of  the  physics  and 
data  processing  techniques  used  in  each  step  and  hence  tends  to  lend 
support  to  the  worth  of  the  final  product  more  so  than  the  performance 
of  only  wave  comparisons,  regardless  of  how  extensive  these  compari¬ 
sons  may  be.  Indeed,  if  each  step  is  shown  to  be  physically  valid,  it 


can  be  argued  that  the  results  should  be  as  accurate  in  sites  where 
there  are  no  wave  data  for  verification  as  they  are  in  areas  where 
large  amounts  of  gage  data  are  available.  Additionally,  if  all  steps 
are  modeled  correctly,  factors  such  as  direction  and  angular  spreading, 
which  are  not  generally  available  for  comparisons,  can  reasonably  be 
assumed  to  be  at  least  approximately  correct. 

This  report  will  discuss  the  numerical  wave  model  used  to  gener¬ 
ate  the  deepwater  wave  data  that  are  contained  in  the  Atlantic  Phase  I 
and  Phase  II  data  reports.  A  FORTRAN  listing  is  included  and  a  dis¬ 
cussion  of  how  to  use  the  model  is  provided  herein. 
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A  NUMERICAL  MODEL  FOR  WIND-WAVE  PREDICTION  IN  DEEP  WATER 


Introduction 

1.  This  report  provides  a  listing  and  description  of  a  simpli¬ 
fied  version  of  the  numerical  wave  hindcast  model  used  by  the  U.  S.  Army 
Engineer  Waterways  Experiment  Station  to  develop  wave  climates  for  U.  S. 
coastal  regions.  The  report  will  outline  the  program  architecture,  time 
and  storage  requirements,  and  input  and  output  parameters.  Also  in¬ 
cluded  is  a  sample  setup  of  a  run  with  a  sample  of  the  output.  It  is 
the  intent  of  this  report  to  provide  a  sufficient  description  of  the 
model  to  allow  a  person  reasonably  knowledgeable  in  numerical  modeling 
of  water  waves  to  run  the  model  with  little  difficulty.  Due  to  the 
complexity  of  the  physics  of  the  processes  modeled  and  the  modeling 
techniques  employed,  it  is  not  possible  to  develop  a  code  that  is  com¬ 
pletely  foolproof;  indeed,  it  is  possible  through  apparently  subtle 
misspecification  of  wind  fields  or  other  input  parameters  to  create 
erroneous  results.  A  description  of  the  physics  of  wave  growth  incor¬ 
porated  in  this  model  is  provided  by  Resio  (1981). 

2.  The  model  presented  here  is  similar  to  the  one  used  to  hind- 
cast  20  years  of  wave  data  in  the  Atlantic  Ocean  (Corson  et  al.  1981). 
However,  the  latter  model  was  written  explicitly  to  operate  on  an 
orthogonal  grid  inscribed  on  a  sphere,  and  much  of  the  code  involved 
nonstandard  FORTRAN.  The  present  model  has  been  modified  to  consider 

a  flat  earth  and  an  attempt  has  been  made  to  code  the  program  in  a  man¬ 
ner  that  should  run  with  little  or  no  changes  on  most  computers.  Appen¬ 
dix  B  provides  a  discussion  of  the  FORTRAN  changes  in  the  coding  to 
convert  the  flat  earth  wave  model  into  a  spherical  orthogonal  system  if 
large  oceanic  areas  are  to  be  modeled.  The  wave  model  requires  a  rect¬ 
angular  grid  and  wind  input  at  each  of  the  intersections  of  the  grid. 

3.  The  model  is  relatively  simple  to  use.  However,  the  user  must 
realize  that  this  wave  model  is  but  one  of  a  sequence  of  models  used  to 
develop  the  wave  climate.  This  sequence  consists  of  programs  to  convert 


pressure  fields  into  wind  fields,  wind  fields  into  wave  fields,  and 


wave  data  into  statistical  summaries.  Considerable  care  has  been  taken 
to  assure  that  all  elements  of  the  models  are  compatible  and  the  models 
have  been  verified  separately  and  together. 

4.  The  following  constraints  are  brought  to  the  attention  of  the 


a.  Wind  input.  The  model  requires  the  input  of  a  wind  ve¬ 
locity  (speed  and  direction)  at  each  grid  point  as  a 
function  of  time.  The  wave  growth  terms  in  the  model  are 
effectively  scaled  on  the  local  wind  stress  based  on  the 
assumption  of  a  neutral  (0°  air-sea  temperature  differ¬ 
ence)  atmospheric  stability.  The  wind  model  that  feeds 
wind  data  to  the  wave  model  accounts  for  both  atmospheric 
stability  and  baroclinicity .  The  user  must  therefore 
transform  observed  or  predicted  surface  wind  fields  to 
give  this  effective  wind  velocity.  The  procedure  used  in 
the  Wave  Information  Study  (WIS)  is  provided  in  WIS  Re¬ 
ports  4  and  10  (in  preparation).  The  user  is  specifically 
warned  that  the  input  of  wind  velocities  such  as  those 
read  directly  from  a  weather  map  into  the  model  without 
making  these  adjustments  can  lead  to  spurious  results. 

b.  Grid  geometry.  The  model  was  designed  for  water  bodies 
of  large  and  relatively  uncomplicated  geometry.  The 
model  employs  a  fourth  order  finite-difference  scheme  for 
most  wave  propagation.  In  the  case  of  an  irregular  shore¬ 
line  combined  with  a  narrow  body  of  water,  most  of  the 
computation  points  are  not  fourth  order  points  as  later 
described,  and  the  model  may  become  numerically  unstable. 
This  model  has  only  been  tested  for  a  simple  geometry. 

c.  Input  parameters.  Constraints  on  numbers  of  frequencies, 
size  of  space,  and  time-steps  must  be  adhered  to  to  avoid 
instabilities  or  inadequate  resolution  of  the  energy 
spectrum. 

d.  Program  modifications.  The  numerical  procedures  and  ap¬ 
proximations  to  the  wave  growth  and  decay  processes  are 
complex.  As  a  result  of  a  desire  to  obtain  computational 
efficiency,  certain  parts  of  the  calculations  are  split 
into  various  parts  of  the  program.  Such  programming 
steps  therefore  make  the  code  less  penetrable  to  some¬ 
one  not  intimately  familiar  with  the  program  development. 
Consequently,  great  care  should  be  taken  in  modifying  any 
part  of  the  program  except  the  input  and  output  portions 
of  the  model. 

5.  Again,  it  is  presumed  that  the  user  of  the  model  is  basically 
familiar  with  the  modeling  of  the  growth  and  decay  of  a  directional 
wave  spectrum.  Therefore,  this  report  is  written  with  a  minimum  of 
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explanations  to  this  background  information.  The  user  not  so  familiar 
may  need  to  make  reference  to  WIS  Reports  1-11. 

Description  of  Computer  Program 

6.  This  section  contains  a  description  of  the  computer  program 
and  includes  flow  charts  and  notes  on  the  dimension  statements  needed 
for  the  variables  in  the  computer  program.  General  aspects  of  running 
the  program  are  included  in  this  section.  The  sections  on  data  input 
and  output  are  self-explanatory  and  are  necessary  to  run  the  program. 
Frequency  and  direction 

7.  The  model  uses  20  discrete  frequencies  and  16  direction  cate¬ 
gories,  giving  320  frequency-direction  components  at  each  spatial  point. 
These  frequencies  are  read  in  by  the  input  file.  Twenty  frequencies 
are  needed  to  run  this  program.  The  numerical  values  for  these  frequen¬ 
cies  can  be  changed,  and  each  Af*  (difference  between  frequencies)  is 
calculated  within  the  program.  It  is  possible  with  minor  modifications 
to  change  the  number  of  frequencies  considered;  however,  the  accuracy 

of  the  wave-wave  interaction  source  term  used  in  this  model  depends  on 
the  resolution  of  the  spectral  peak.  If  too  crude  a  representation  is 
used  to  define  the  spectrum,  problems  arise  in  the  estimation  of  these 
source  terms.  The  number  of  direction  categories  cannot  be  changed 
without  modifying  the  whole  program.  The  part  of  the  spectrum  at 
frequencies  higher  than  the  highest  frequency  input  is  calculated 
parametrically. 

Core  storage 

8.  Core  storage  requirements  are  minimized  by  performing  all 
calculations  for  each  frequency  separately.  In  this  way,  only  one 
large  matrix  containing  all  of  the  two-dimensional  spectral  information 
is  needed.  The  core  storage  required  to  load  this  program  is  described 
in  Figure  1  in  terms  of  the  number  of  spatial  points  (water  points, 


*  For  convenience,  symbols  and  unusual  abbreviations  are  listed  and 
defined  in  the  Notation  (Appendix  D). 
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GRID  SIZE  RELATED  ARRAYS 


VARIABLE  DIMENSIONS 


»  M  *  ROWS  x  COLUMNS  IN  GRID 
S  •  NUMBER  OF  GRID  POINTS 
IN  »  M  -  SI 

LW  -  NUMBER  OF  WATER  POINTS 
BB  •  NUMBER  OF  BOUNOARY  POINTS 
ILW  *  BB  I  SI 


FSWLA(S),  SWLSC(S) 
HSCALE(S),  ESUM(S),  VSUM(S) 
USUM(S),  ALPHA(S),  EMXA(S), 
FOAIS).  THOA(S),  ISWAR(S), 
EMAXSW(S),  FSUMIS),  EF(S), 
SCAA(S),  FOOOIS),  EA(S) 
COMMON/EE/  LOCKS),  LOCJIS) 
COMMON/NP/  KFAfSI,  MLP(S) 
IASWIS),  ESWSQ(S),  ISWO(S), 
AF07IS) 


IAMT(S,16),  IBMTIS.16) 
ICMT(S,16) 

COMMON  /RP/  KSSULW(LW,  16) 
COMMON  /KU/  KSSUA  (BB,  16) 
COMMON  /A/  E  (S,16),  ENA(S,16) 


I  INGLEI16),  COSRI16) 

|  SINRI16),  GAMANGI361) 

|  EPHON(16),  IAAD(361) 
COMMON  /IJ/  IUAI16),  JUAI16) 


DF2I20),  FF(21),  ANGL(20) 

C( 20) ,  XCRI21),  DELF(20) 

SINF(20),  RM33I21) 

RM22I21),  IAADW(20) 

COMMON  /DD/  CG(20),  DELX,  TINC 
COMMON  /ST/  PUK20),  PU2I20) 
F11(21),  ZFLI21),  FSEAI20.21) 
EFSUMI20),  KREFA(21) 


LOUT  (20,2) 

COMMON  /B /  PMU  (5.20) 


TEL  (S,  16,20) 


I  COMMON  /ST/  XMA1(  16,20), 

I  XMA2  (16,20) _ 

3!  COMMON  /8/  VMU  (5,16,20) 


COMMON  /RP/  EVI5.16),  LXLK16.5), 
LXU(16,5),  IREFAI16) 
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WINDA(N.M),  WDIR(N,M),  H2(N,M) 
COMMON  /CB/  LA(N,M),  LOCC(N,M 

I - * 

|  EPMTRX  (5,S)  6 

!"cOSSQ(181l.  SNLA(40C) 

1  FMX(IOO),  FFLIMdOOl 
I  USTA(IOO),  SNLSW(400) 

'J_TITL(2),  AFQ7(100)  13 


MISCELLANEOUS  DIMENSIONS  THAT  DO  NOT 
CHANGE  WITH  SIZE  OF  GRID 


Figure  1.  Description  of  core  storage  and  dimension  values  needed  fo 
wave  model  program  arrays.  Variable  dimensions  (S,  N,  M,  LW,  BB)  can 
be  easily  changed  for  different  grid  sizes,  but  a  change  of  the  numer 
leal  dimensions  alters  the  program.  The  arrays  are  grouped  by  types 

of  arrays 


DIRECTION  RELATED  VARIABLES 


boundary  points,  and  total  grid  points),  the  number  of  frequencies,  the 
number  of  direction  rays,  and  the  wind  speeds  used  in  the  calculations. 
The  grid  size  related  arrays  are  the  only  ones  that  need  to  be  modified 
for  a  new  size  grid.  The  frequency,  wind  speed,  and  direction  related 
variables  stay  constant. 

9.  To  set  up  the  dimension  statements  for  a  specific  grid,  first 
determine  the  size  of  the  grid  using 
N  =  number  of  rows  in  grid 
M  =  number  of  columns  in  grid 

where  N  and  M  are  chosen  to  portray  an  adequate  resolution  of  the 
water  body  modeled.  (N  and  M  are  numbered  in  the  row-column  conven¬ 
tion  of  a  matrix  where  (1,  1)  is  in  the  upper  left-hand  corner.)  The 
number  of  grid  points  (S)  is  equal  to  N  times  M  .  Second,  since  this 
grid  is  a  map  representation  of  some  body  of  water,  it  is  necessary  to 
determine  which  of  the  grid  point  intersections  are  land  points  and 
which  are  water  points.  Some  judgment  is  needed  to  set  up  this  grid 
to  portray  a  smooth  body  of  water.  Small  islands  or  points  of  land 
jutting  out  into  the  water  may  have  to  be  assumed  to  be  water  points 
even  though  a  grid  intersection  may  actually  have  land  under  it.  A 
sample  input  grid  is  shown  in  Appendix  C.  Third,  determine  the  number 
of  Lax-Wendroff  water  points  (LW)  and  the  number  of  boundary  points 
(BB).  A  true  Lax-Wendroff  water  point  (fourth-order  point)  is  defined 
as  having  two  water  points  surrounding  it  in  every  direction.  The 
boundary  points  will  include  the  other  points  that  are  not  Lax-Wendroff 
points.  LW  plus  BB  should  be  less  than  or  equal  to  S  .  It  is 
possible  to  dimension  the  model  using  only  water  points  rather  than 
the  total  number  of  grid  points  for  the  various  source  terms.  This 
will  save  on  the  amount  of  core  storage  if  the  grid  is  very  large. 

This  report  describes  the  method  of  using  the  S  variable  as  the  total 
size  of  the  grid. 

10.  The  total  amount  of  core  storage  can  be  computed  using  S  , 

N  ,  M  ,  LW  ,  and  BB  and  the  equation  in  Figure  1.  For  most  grids, 
it  is  easier  to  let  LW  =  BB  =  S  where  S  =  N  x  M  .  For  large  grids, 
let  S  =  LW  =  BB  =  number  of  water  points.  The  amount  of  core  storage 


is  proportional  to  the  grid  size  (S) .  Figure  1  will  be  discussed  again 
in  the  section  on  model  structure. 

Run  time 

11.  One  of  the  problems  with  writing  a  generalized  code  is  that 
the  program  cannot  be  designed  to  be  optimal  in  terms  of  run  time  on  a 
particular  computer.  Figure  2  provides  a  set  of  comparative  run  times 
on  two  different  computers  for  the  test  run  described  in  Appendix  C. 

Central  Processor  Time*  _ Name  of  Computer _ 

6.2  sec  Cray  Computer**  at  Boeing 

Computer  Services 

348.0  sec  Honeywell  DPS-1  at  WES 

*  The  computer  times  represented  here  refer  to  the  run  in 
Appendix  C  which  has  a  9*9  grid  and  a  grid  spacing  of 
222.2  km.  The  time-step  was  3,600  sec.  Seventy-two  time- 
steps  were  taken. 

**  Note  that  the  computer  listed  here  is  a  vectorizing 
machine.  The  code  has  been  designed  for  a  vectorizing 
machine  and  run  times  will  increase  significantly  if 
vectorizing  does  not  take  place. 

Figure  2.  Comparative  run  time  on  different  computers 
Model  structure 

12.  Figures  3-5  show  the  basic  structure  of  the  model.  After  an 
initialization  section  of  code,  the  calculations  are  completely  con¬ 
tained  within  a  series  of  loops.  Besides  setting  constants  and  precal¬ 
culating  various  spectral  shape  functions  and  directional  distributions, 
the  initialization  section  of  code  calls  the  subroutine  CLASSB  which  de¬ 
termines  the  referencing  counters  for  all  spatial  points  along  with  the 
propagation  area  (Lax-Wendrof f  or  upstream  differencing)  in  which  they 
lie.  CLASSB  also  calculates  the  propagation  multipliers  to  be  used  in 
subroutine  PROPR,  which  performs  the  actual  propagation  calculations. 

13.  Appendix  A  is  a  FORTRAN  listing  of  the  wave  model  program. 
Comments  have  been  added  to  the  program  to  give  the  user  an  idea  what 
section  of  code  corresponds  to  what  physical  process.  The  example  in 
Appendix  C  was  run  on  a  Cray  computer  system.  If  the  grid  dimensions 
are  different  from  the  example  model  and  core  storage  is  minimal, 
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INPUT  SETS  1  -  4 


INITIALIZATION  OF  CONSTANTS  AND  PRECALCULATIONS 
OF  PARAMETERS  INDEPENDENT  OF  TIME 


CALL  CLASS6  (FIGURE  4) 

- - OPTIONAL  RESTART  INFORMATION  READ 

* - OPTIONAL  SETUP  OF  TEST  WIND  CONDITIONS 

- —  _ 

TIME  LOOP 

I - OPTIONAL  WIND  INPUT 


ESTIMATION  OF  NONFREQUENCY  DEPENDENT  PARAMETERS 


OUTPUT  ROUTINES 


OPTIONAL  OUTPUT  OF  RESTART  PARAMETERS 

END 


Figure  3.  Fundamental  structure  of  program 


DETERMINE  MAPPING  INTO  WATER-ONLY  REFERENCING  SYSTEM  SEPARATING 
BOUNDARY  POINTS  AND  INTERIOR  POINTS  INTO  DIFFERENT  SEQUENCES 

♦ 

CALCULATE  PROPAGATION  MULTIPLIERS  FOR  UPSTREAM  DIFFERENCING  REGiON 

I 

CALCULATE  PROPAGATION  MULTIPLIERS  FOR  LAX-WENDROFF  REGION 

* 

DETERMINE  LOCATION  REFERENCES  FOR  UPSTREAM  DIFFERENCE  POINTS 

I 

DETERMINE  LOCATION  REFERENCES  FOR  INTERIOR  POINTS 


Figure  4.  Structure  of  subroutine  CLASSB 


PROPAGATE  ENERGIES  AT  INTERIOR  POINTS  (PROPAGATION  METHOD  CAN 
BE  EITHER  LAX-WENDROFF  OR  UPSTREAM  DIFFERENCING  DEPENDING  ON 
LOCATION  OF  FREQUENCY  RELATIVE  TO  LOCAL  SEA  PEAK) 


-  COMPLETE  PROPAGATION  ALONG  8  PRIMARY  AXES 


-  FIRST  PART  OF  PROPAGATION  ALONG  8  OFF-AXIS  DIRECTIONS 

♦ 

PROPAGATE  ENERGIES  AT  BOUNDARY  POINTS 


-  COMPLETE  PROPAGATION  ALONG  8  PRIMARY  AXES 


-  FIRST  PART  OF  PROPAGATION  ALONG  8  OFF-AXIS  DIRECTIONS 

I 

PROPAGATE  ENERGIES  AT  INTERIOR  POINTS  TO  COMPLETE  PROPAGATION 
FOR  OFF-AXIS  DIRECTIONS  I 


PROPAGATE  ENERGIES  AT  BOUNDARY  POINTS  TO  COMPLETE  PROPAGATION 
FOR  OFF-AXIS  DIRECTIONS  | 


CHECK  ON  NEGATIVE  ENERGY  AT  EACH  POINT  (E  (f  ,6)  >  0) 


Figure  5.  Structure  of  subroutine  PROPR 
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changes  need  to  be  made  in  the  variables  S  ,  N  ,  M  ,  LW  ,  and  BB 
(see  Figure  1)  to  match  the  specific  grid  being  considered.  The  FOR¬ 
TRAN  listing  in  Appendix  A  is  dimensioned  large  enough  to  accept  a  31x31 
grid  and  requires  about  520K  words  of  storage.  The  dimensions  also  need 
to  be  changed  in  the  two  subroutines  (CLASSB  and  PROPR)  to  match  the 
main  program's  dimensions. 

14.  The  model  can  begin  calculations  from  a  cold-start  condition 
(i.e.,  initialize  everything  to  preset  values)  or  it  can  begin  from  a 
warm-start  condition  with  initial  spectral  and  parametric  values  being 
read  in  from  a  file.  An  option  to  create  such  a  restart  file  is  in¬ 
cluded  at  the  end  of  each  run.  If  this  code  is  used  to  calculate  a 
long  time  series,  it  is  advisable  to  break  the  series  into  smaller  in¬ 
crements  and  run  these  increments  sequentially  using  the  restart  capa¬ 
bility  of  the  model. 


Input  to  the  Model 


15.  The  input  and  output  for  this  code  were  kept  simple  to  main¬ 
tain  generality  and  computer  compatibility.  The  wave  model  input  con¬ 
sists  of  a  set  of  data  input  cards  that  describe  the  grid,  options, 
frequencies,  and  initial  conditions  and  one  set  of  data  that  describes 
the  wind  in  terms  of  time,  speed,  and  direction.  Program  modifications 
to  produce  special  input  or  output  routines  can  be  made  by  users  with 
specific  requirements.  However,  such  routines,  especially  computer 
graphics,  tend  to  vary  from  computer  to  computer.  Therefore,  such 
options  have  been  left  to  the  discretion  of  the  user. 

Data  input 

16.  Data  input  for  this  model  is  contained  on  four  sets  of  data 
cards  (these  four  card  sets  were  combined  into  input  file  INPTEST  in 
Appendix  C) : 

a .  Card  Set  1:  Problem  Description  (1  card) 

Input  Parameters  _ Description _ 


N 

M 


Number  of  rows  in  grid 
Number  of  columns  in  grid 


TINC*  Number  of  seconds  per  time-step 

DINC  Distance  between  grid  points  in 

kilometres  (DELX) 

MSTA  Number  of  special  output  locations 

IHR  Number  of  hours  between  wind  inputs 

NEWSTR  1  for  cold  start;  0  for  warm  start 

(NEWSTR  should  be  set  to  1  for  test 
runs  or  for  runs  without  restart 
data  from  previous  run) 

NOWRT  1  for  no  tape  writes;  0  for  output 

tape  writes  (NOWRT  should  be  set 
to  1  for  test  runs) 

NOR D  1  for  no  tape  reads;  0  for  input 

tape  reads  (NORD  should  be  set  to 
1  for  test  runs) 

ZM2  Anemometer  level  of  wind  input  in 

centimetres 


*  The  time-step  is  constrained  by  the  stability  criterion 

for  propagation  of  the  lowest  frequency,  i.e.,  TINC  <  DELX/ 

c  (f  .  )  .  DELX  (the  X-value  of  DINC)  is  distance  between 
g  min 

grid  points  and  Cg(^min)  *s  t*le  8rouP  velocity  of  the  mini¬ 
mum  frequency  band  that  is  input  into  the  wave  model. 

b.  Card  Set  2;  Discretized  Frequencies  (2  cards) 

These  two  data  cards  contain  the  20  frequencies  selected 
for  representing  the  wave  spectrum. 

c.  Card  Set  3:  Special  Output  Locations  (MSTA  cards) 

These  MSTA  pairs  of  numbers  (one  pair  per  card)  are 
the  I,  J,  coordinates  of  the  special  output  locations 
desired.  It  should  be  noted  here  that  if  MSTA  equals 
zero,  no  data  card(s)  will  be  read. 

d.  Card  Set  4:  Land-Sea  Grid  (N  cards) 

This  final  set  of  N  data  cards  contains  the  land  or  water 
data  (input  by  row)  for  the  grid  being  used.  A  value 
<  0  represents  land  and  a  value  >  0  denotes  water.  (This 
input  matrix  of  land-water  points  should  contain  a  border 
of  1  row  and  1  column  of  land  points.) 


17.  Most  of  the  above  variables  are  self-explanatory.  The  num¬ 
ber  of  special  output  locations  is  the  number  of  grid  intersections  that 
are  of  special  interest  to  the  user.  The  wave  information  for  these 


f 

points  will  be  printed  out  and  labeled  separately. 

18.  The  anemometer  level  of  the  wind  input  for  this  wave  model  is 
read  in  as  an  input  parameter,  with  friction  velocities  used  to  scale 
various  parameters  inside  the  program.  A  misspecification  of  this  input 
level  could  lead  to  significant  biases  in  wave  predictions.  For  exam¬ 
ple,  if  10-m  winds  are  assumed  and  the  actual  level  of  the  wind  field 
was  taken  at  19.5  m,  calculated  waves  would  be  consistently  too  high. 

Wind  input 

19.  As  discussed  in  the  introduction,  wind  speeds  coming  into 
the  wave  model  are  transformed  using  the  following  theoretical  consid¬ 
erations.  The  wind  profile  equation,  used  to  obtain  friction  velocity 
from  the  wind  velocity  inside  the  wave  model,  assumes  a  particular  drag 
law  as  well  as  neutral  stability.  The  drag  law  used  is  consistent  with 
that  obtained  by  Garratt  (1977)  for  velocities  over  12  m/sec.  Liu  and 
Ross  (1980)  have  shown  that  stability  has  a  strong  influence  on  wave 
growth.  Consequently,  wind  speeds  coming  into  the  wave  model  are  taken 
to  be  equivalent  neutral  wind  velocities,  i.e.,  those  wind  velocities 
at  neutral  stability  which  produce  the  same  friction  velocity  as  the 
actual  recorded  or  estimated  winds  in  a  given  nonneutral  condition. 

The  methods  used  by  WIS  are  given  in  WIS  Reports  4  and  10  (in 
preparation) . 

20.  The  form  of  the  wind  input  for  the  program  will  be  left  to 
the  discretion  of  the  user.  The  input  file  ILWND  in  Appendix  C  is  an 
example  of  a  wind  input  file.  The  user  is  reminded  that  wind  input  is 
the  primary  driver  of  the  wave  model.  The  WIS  wind  data  came  from  a 
wind  model  that  took  existing  pressure  data  and  converted  this  data  to 
a  wind  field  (Corson,  Resio,  Vincent  1980;  Resio,  Vincent,  Corson 
1982).  Any  extreme  variations  in  the  wind  field  will  be  reflected  in 
unusual  wave  conditions.  If  unusual  wave  conditions  appear  as  wave 
model  output,  the  user  is  advised  to  check  the  input  wind  field  first. 

The  program  contains  an  option  (NORD  =  1)  that  allows  the  user  to  enter 
a  constant  wind  speed  and  direction  for  test  purposes.  These  constants 
can  be  entered  on  a  card  immediately  following  card  set  4  in  the  format 


specified  in  the  program.  It  is  suggested  that  the  user  try  a  test  run 


using  constant  winds  before  attempting  a  run  using  an  input  wind  tape. 
Wind  speeds  are  entered  in  knots  and  wind  direction  is  in  degrees.  The 
wave  model  will  not  accept  zero  wind  speeds  but  a  zero  in  the  direction 
category  is  acceptable.  Wind  data  should  be  in  integer  format.  Using 
the  NORD  =  0  option,  a  wind  data  tape  can  be  used  as  input.  The  wind 
data  tape  needs  to  have  a  specification  of  date  and  hour  in  integer 
form.  The  wind  speeds  and  directions  are  read  in  for  each  point  of  the 
grid  at  the  specified  time  intervals.  The  directions  of  the  input 
winds  should  be  referenced  to  the  direction  toward  which  they  are  blow¬ 
ing.  The  direction  system  in  the  wave  model  is  a  polar  co-ordinate 
system  with  0°  corresponding  to  east  and  90°  corresponding  to  north. 

A  0°  direction  for  a  wind  indicates  a  wind  blowing  from  west  to  east. 
Frequency  input 

21.  The  discretized  frequencies  can  be  varied  by  the  user.  Sug¬ 
gested  values  in  Hz  for  the  20  frequencies  for  oceanic  conditions  are: 
0.03,  0.04,  0.05,  0.06,  0.07,  0.08,  0.09,  0.10,  0.11,  0.12,  0.13,  0.14, 
0.15,  0.16,  0.17,  0.18,  0.19,  0.20,  0.21,  and  0.22.  The  time-step 

is  constrained  by  the  group  velocity  of  the  minimum  frequency 
[TINC  <  DELX/cg(fmin)]  ;  therefore,  in  situations  where  the  grid  spac¬ 
ing,  DELX,  is  small  and  the  time-steps  are  very  small,  fmin  can 
increased  to  create  larger  time-steps  (less  computer  time).  The 
difference  between  frequencies,  Af  ,  can  remain  constant  as  shown  or 
vary  throughout  the  input  frequency  range.  The  frequency  range  must 
be  realistic  for  the  spectral  situations  being  considered. 

Standard  Output  from  the  Model 

22.  The  output  of  the  model  is  divided  into  three  types.  (See 
Appendix  C  for  an  example  of  the  output.)  The  first  part  of  the  output 
is  a  listing  of  the  input  parameters  such  as  the  type  of  start  (cold  or 
warm),  the  grid  size  (NxM),  the  time  increment  (TINC)  in  seconds,  and 
the  grid  spacing  (DINC)  in  kilometres,  the  input  frequencies,  and  the 

I  and  J  grid  location  of  the  special  output  stations.  The  grid 
geometry  shows  the  actual  grid  that  was  input  and  it  also  shows  a  grid 


with  the  KSS  number  for  each  of  the  water  points  of  the  grid.  The  inte¬ 
rior  numbers  beginning  with  one  show  the  Lax-Wendroff  (fourth  order) 
points.  The  second  group  of  numbers  on  this  grid  count  the  boundary 
water  points  that  will  be  used  in  upstream  differencing.  The  numbers 
below  the  grid  refer  to  the  total  number  of  water  points,  the  number  of 
fourth-order  points,  and  the  number  of  boundary  points.  If  a  constant 
wind  test  was  run  for  test  purposes,  the  test  parameters  appear  below 
the  grids  showing  the  input  wind  speed,  wind  direction,  and  the  number 
of  time-steps.  The  second  type  is  the  information  needed  at  the  spe¬ 
cial  output  locations  defined  in  the  input.  This  output  is  discussed 
in  paragraph  23.  The  third  type  is  the  information  at  each  point  of 
the  entire  grid.  This  output  is  discussed  in  paragraph  24.  The  infor¬ 
mation  contained  in  the  special  output  prints  can  be  obtained  for  every 
point  of  the  grid;  but  if  only  a  few  sites  are  needed,  it  is  necessary 
only  to  print  this  out  for  a  few  locations.  Calculations  are  done 
using  centimetre-gram-second  units  and  numerical  values  are  in  these 
units  unless  otherwise  noted. 

23.  The  output  for  the  special  output  locations  is  printed  for 
each  time-step.  The  date-time  and  time  in  hours  since  start  of  run 
precedes  this  output.  Each  special  output  location  is  labeled  by  its 
I  and  J  location  (matrix  row-column  convention)  in  relation  to  the  grid. 
The  first  type  of  information  included  for  the  special  output  location 
includes  a  listing  of  the  spectral  energy  by  frequency  and  angle  loca¬ 
tion.  The  one-dimensional  energy  integrated  over  the  16  directions  is 
given  per  frequency  in  the  column  headed  by  Energy  Density.  The  rest 
of  the  chart  lists  the  energy  by  frequency  in  the  direction  band  in 
which  it  occurs.  The  special  output  information  below  the  chart  is 
defined  by  its  title  and  units.  The  total  energy  is  the  energy  in  the 
discrete  frequency  range  and  the  parametric  energy  is  the  energy  in  the 
frequencies  above  the  discrete  region.  The  significant  wave  height  is 
calculated  by  multiplying  0.0401  times  the  square  root  of  the  total 
energy.  The  rest  of  the  values  are  self-explanatory.  By  listing  no 
special  output  stations,  all  the  special  output  printouts  can  be 
deleted. 
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24.  The  output  for  the  entire  spatial  grid  is  printed  in  matrix 
fora.  The  set  of  matrices  is  preceded  by  a  date- time  heading.  The 
first  matrix  contains  the  significant  wave  height.  The  significant 
wave  height  in  metres  is  displayed  for  all  grid  points  at  each  time- 
step  and  is  calculated  by  multiplying  0.0401  times  the  square  root  of 
the  total  energy.  The  data  can  be  located  by  the  I  and  J  locations  on 
the  spatial  grid.  The  next  matrix  contains  the  mean  direction  of  the 
local  sea  and  the  third  matrix  contains  the  period  of  the  local  sea. 
The  printout  of  significant  wave  height  presented  in  Appendix  C  is 
typical  of  these  general  output  matrices.  Output  tape  writes  can  be 
added  if  needed  by  the  user. 

Model  Characteristics 

25.  An  exaiaple  of  a  complete  wave  model  run  for  a  variable  wind 
field  is  given  in  Appendix  C.  Although  such  a  run  serves  to  demon¬ 
strate  a  number  of  different  aspects  of  the  wave  generation,  propaga¬ 
tion,  and  decay  characteristics  of  the  wave  model,  it  does  not  provide 
a  clear  understanding  of  the  overall  model  characteristics.  To  accom¬ 
plish  this,  two  types  of  idealized  generation  cases  are  presented  in 
this  section — fetch-limited  wave  growth  and  duration-limited  wave 
growth  up  to  a  fully  developed  sea. 

26.  Figure  6  shows  a  comparison  of  the  rates  of  growth  produced 
by  the  wave  model  and  the  theoretical  rates  of  growth  from  Hasselmann 
(1976).  The  wave  model  reproduces  the  observed  fetch  limited  growth 
patterns  quite  well.  Figure  7  gives  the  rate  of  growth  of  waves  as  a 
function  of  time.  It  is  apparent  from  this  figure  that  the  rate  of 
wave  growth  begins  to  be  significantly  altered  at  a  nondimensional 
time  (t  =  gt/u)  of  about  2.1  x  106.  The  approach  to  fully  developed 
conditions  is  not  accomplished  by  constraining  the  wave  spectrum  to  a 
predetermined  fully  developed  spectral  form,  but  is  done  by  limiting 
the  scaling  frequency  for  wave-wave  interactions  to  a  value  dependent 
on  the  wind  speed.  Thus,  a  steady-state  condition  is  never  actually 
reached;  instead,  the  rate  of  growth  asymptotically  approaches  zero. 
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27.  Three  tests  should  be  run  to  verify  the  boundary  conditions 
and  stability  of  the  wave  model.  These  tests  consist  of  a  fetch  test 
(where  the  results  can  be  compared  with  Figure  6),  a  time  growth  test 
(where  results  can  be  compared  with  Figure  7),  and  a  test  where  the  wind 
field  shifts  by  90  deg.  Each  test  should  be  run  for  several  different 
constant  winds  such  as  20,  30,  40,  and  60  knots.  The  fetch-limited  test 
should  be  run  for  30  hours,  and  the  duration  test  should  be  run  for 

150  hours.  The  wave  model  has  been  set  up  with  an  option  to  run  with 
a  file  containing  constant  winds  (where  each  grid  point  has  the  same 
value  of  wind  speed  and  direction).  This  corresponds  to  using  a  value 
of  one  for  the  NORD  variable  described  in  paragraph  16.  The  user  can 
add  the  necessary  logic  to  the  program  to  create  the  fetch-limited  and 
the  wind  shift  conditions.  These  tests  have  been  run  using  the  FORTRAN 
in  Appendix  A.  If  any  significant  changes  are  made  in  the  program, 
these  tests  need  to  be  run  again. 

Discussion 

28.  The  wave  model  described  in  this  report  is  a  state-of-the-art 
method  of  using  wind  data  to  hindcast  wave  heights  in  deep  water.  The 
model  uses  the  equations  that  depict  the  physics  of  the  wave  interac¬ 
tion  and  growth  and  uses  numerical  procedures  to  evaluate  the  results 

of  these  equations.  This  model  has  been  termed  a  discrete  spectral 
model  (Resio  1981)  and  models  the  spectrum  by  using  a  discrete  band  of 
frequency-direction  components  on  the  forward  face  of  the  spectrum. 

Other  reports  in  the  Wave  Information  Study  deal  with  the  confidence 
levels  of  the  numerical  results  and  compare  calculated  results  with 
actual  gage  data  (Corson  and  Resio  1981,  Baird  and  Readshow  1981). 

The  program  has  been  written  in  a  form  that  will  be  compatible  on 
several  different  computer  systems  and  is  especially  appropriate  for 
machines  that  have  a  large  core  storage  capacity  and  the  potential 
for  vectorization. 

29.  As  stated  in  the  introduction,  this  code  does  not  contain 
the  changes  needed  for  a  spherical  orthogonal  grid  system.  This  code 
has  been  used  in  areas  such  as  bays  and  harbors  and  small  sections  of 


the  ocean.  Deep  water  is  assumed  throughout  all  calculations. 

30.  The  validity  of  the  wave  data  generated  with  the  WIS  wave 
model  is  directly  related  to  the  accuracy  of  the  input  wind  data.  For 
WIS  hindcasts,  every  effort  is  made  to  generate  representative  wind 
fields.  The  procedures  used  to  generate  accurate  wind  fields  and  com¬ 
parisons  of  hindcast  and  measured  winds  are  discussed  in  Resio,  Vincent, 
Corson  (1982).  In  order  to  assure  valid  wave  computations,  users  of 
the  WIS  wave  model  should  verify  input  winds  thoroughly. 
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APPENDIX  A:  WAVE  MODEL  LISTING 


The  following  is  a  FORTRAN  listing  for  the  flat  bed  wave  model. 
The  "warm  start"  feature  has  been  included  in  the  proper  place  but  has 
been  entered  as  comment  statements  to  allow  the  user  to  modify  this  for 
the  computer  system  he  will  be  using.  This  FORTRAN  listing  was  used  to 
produce  the  sample  output  included  in  Appendix  C  of  this  report. 


noon  oooooorjoo 


PROGRAM  WMRWAVE 

C  UMRUAVE  IS  A  PHASE  I  FLAT  BED  HAVE  NOBEL 

C  THIS  VERSION  HAS  DESIGNED  FOR  THE  HAVE  MODEL  REPORT  (1982) 

C  DIMENSIONS  HILL  ACCOMODATE  961  HATER  POINTS  ON 
C  A  31X31  GRID. 

DIMENSION  FSULA(961 ) » SULSC (961 ) 

DIMENSION  IPP<961) 

DIMENSION  INGLEU6) 

DIMENSION  LOUT ( 20 . 2 ) . DF2 ( 20 ) . FF< 21 ) . ANGL ( 20 )  r 
1C0SR( 16) *SINR< 16) *6AMANG(361 ).C(20).XCR(21 ) * 

2HINDA(31 »31 ) .UDIR(31>31) iHSCALE(961 ) »DELF(20)  >SINF(20) 

DIMENSION  C0SSQ( 181 > >EPMTRX(S> 961 ) >SNLA(400) >RM33(21 ) . RM22(21 ) » 
1FNX(100>.FFLIM<100>.USTA<100).ESUM(961).VSUN(961).USUM(961). 
2EPH0N( 16) > ALPHAC961 ) »EMXA(961 ). H2 (36.36). F0A( 961).  TH0A( 961 ) > 

31 A AD ( 361 ) . I AADW(20) > ISWARC961 ) iEMAXSU<961 ).FSUN<961).EF(961) 
DIMENSION  SCAA(961) iF000(961) 

DIMENSION  I AMT (961.16). IBMT (961 *  16) . ICHT (961 .16) 

COMMON  /DD/CG(20) iDELXiTINC 

COMMON  /ST/XMA1 ( 16i 20) >XMA2< 16. 20) *PU1 (20) >PU2(20) 

COMMON  /CC/  N.M.NPLH.NPNLW.NPTS.NFREQ.NDIR.NPTSPl 
COMMON  /EE/LOCI (961 ) >L0CJ(961 ) 

COMMON  /RP/KSSULU(729i 16) iEV(Si 16) tLXLI ( 1 6 r 5 > *LXLJ( 16.5) 1 1 REF A ( 16) 
C0MM9«  /KU/KSSUA(232* 16) 

COMMON  /IJ/  IUA( 16) .  JUA( 16) 

COMMON  /NP/KFA(961),MLP<961) 

COMMON  /CB/LA(31f31) >L0CC(36i36) 

CGMMON  /B/  VMU (5. 16.20) *PMU( 5.20) 

DIMENSION  EL (961 .20) . TEL (961 . 16.20) 

COMMON  /A/  E( 961  *  16) > ENA <961 *16) 

DIMENSION  F11(21>*ZFL(21) >SNLSH<400) *FSEA(20*21 > 

DIMENSION  EA(961).TITL(2) 

'•IMENSION  IASH(96U  *ESUSG(961 )  *  ISU0(961 ) 

DIMENSION  EFSUN(20) 

DIMENSION  AF07<961 ) . AFQ7< 100) »KREFA(21 ) 

INTEGER  UDIRiUINDA* AF07*AF07*HINDY 
EQUIVALENCE  (NPTS.NPTOT) 

EQUIVALENCE  (TINCtDELT) 

C  PARAMETER  PI=3. 14159*NFREQ=20*NDIR=16*G=980. * TH0PI=6. 28318 
C  PARAMETER  C1=0. 1525. C2=l . A7E-5*C3=0. 00371 *ZM2=1950. *GSQ=960400 
DATA  TITL/4HC0LD.4HUARM/ 

DATA  PI .NFREQ.NDIR.G.THOPI/3. 14159.20. 16*980. *6. 28318/ 

DATA  Cl .C2.C3.GSQ/0. 1525. 1.47E-5.0. 00371. 960400./ 

DATA  IUA/240. 5*1 .3»0. 5*-l .0/. JUA/39-1 . 3*0. 5*1 .3*0.2*-l/ 

DATA  IREFA/O. 5. 0.1. 0.9. 0.5. 0.13. 0.9. 0.1. 0.13/ 

C 

N  =  NUMBER  OF  ROUS  IN  GRID 

M  =  NUMBER  OF  COLUMNS  IN  GRID 

TINC  =  NUMBER  OF  SECONDS  PER  TIME  STEP 

DlNC  =  DISTANCE  BETUEEN  GRID  POINTS  IN  KILOMETERS 

MSTA  =  NUMBER  OF  SPECIAL  OUTPUT  LOCATIONS 
IHR  =  NUMBER  OF  HOURS  BETUEEN  WIND  INPUTS 
NEUSTR  =  COLD  OR  HARM  START 

1  FOR  COLD  START  (NO  RE-START  DATA  FROM  PREVIOUS  RUN) 

0  FOR  WARM  START 

NOURT  =  DETERMINES  WHETHER  OUTPUT  TAPES  ARE  HRITTEN  OR  NOT 
1  FOR  NO  TAPE  WRITES 
0  FOR  TAPE  WRITES 

NORD  =  DETERMINES  WHETHER  INPUT  TAPES  ARE  READ  OR  NOT 


C  1  FOR  NO  TARE  READS 

C  0  FOR  TAPE  READS 

C  ZM2  =  LEVEL  OF  HIND  INPUT  ABOUT  HATER  SURFACE  IN  CENTIMETERS 
C 

NFRQ  =  NFREQ 
NFROl  =  NFREQ  -  1 

C  INITIALIZE  FSULA  FOR  NGR-POINT  GRID 

NGR=  N*H 
DO  713  I=1»NGR 
713  FSWLA( I )=1 .0 

C 

C  READ  PROBLEM  DESCRIPTION  PARAMETERS 

READ  5O4O,N»M,TINC,DINC»MSTA,IHR»NEHSTR»NOHRT»N0RD»ZM2 
5040  FORMAT  (2I5.F10.0.F5. 1 » 5I5»F10.0) 

IF  (NEUSTR.NE.l)  NEUSTR  =  2 
C 

C  PRINT  TYPE  OF  START  (COLD  OR  HARM)  FOR  THIS  RUN 
PRINT  9090.  TITL(NEHSTR) .TITL(NEHSTR) .TITL(NEHSTR) 

9090  FORMAT ( 10X.3(A4. 20X) ) 

C 

C  PRINT  SELECTEB  INPUT  PARAMETERS 

PRINT  2000.  N. M. TINC. DINC 
2000  FORHAT  (/.5X.4HN  =  .I3.3X.4HM  -  .I3.3X. 

.  7HTINC  =  . F8. 0 .3X i 7HDINC  =  .F8.0) 

NDIR2=NDIR/2 

DELX~100000.tDINC 

HR=IHR*3400. 

RTIMES=HR/TINC 

IF(KTIMES*TINC-HR.GT. 0.0001)  GOTO  9999 
C 

C  READ  DISCRETIZED  FREQUENCIES  THAT  REPRESENT  HAVE  SPECTRUM 
READ  220. (FF( I ) . 1=1 i 10) 

READ  220. (FF(I). 1=11.20) 

220  FORMAT  (10F5.3) 

FF(NFREQ+1)=1. 

C 

C  PRINT  DISCRETIZED  FREQUENCIES 
PRINT  800 
800  FORMAT  (/) 

PRINT  8005.  (FF(I). 1=1, NFREQ). 

8005  FORMAT <17H  +++++FREQUENCIES.10F10.5) 

IF(MSTA.LT.l)  GOTO  4450 
C 

C  READ  GRID  COORDINATES  FOR  SPECIAL  OUTPUT  LOCATIONS 
PRINT  800 

DO  118  ISTA=1 »MSTA 

READ  101,  LOUT ( ISTA, 1 ) .LOUT ( I  ST A, 2) 

101  FORMAT  <  4012) 

118  CONTINUE 
C 

C  PRINT  SPECIAL  OUTPUT  LOCATIONS 

PRINT  8009,  ( (LOUT (I,J)»J=1,2),I=1»MSTA> 

8009  FORMAT (20H  SP  OUTPUT  LOCATIONS, 2110) 

4450  CONTINUE 
C 

PRINT  800 

C  READ  AND  PRINT  LAND-SEA  GRID 

C  VALUE  .GT.  0  =  HATER 

C  VALUE  .L£.  0  *  LAND 


C  BATHYMETRIC  DATA  INPUT  BY  ROW 

PRINT  8001 

8001  FORMAT (20H  ++++♦  GRID  GEOMETRY) 

DO  1  1=1 » N 

READ  101 f  <LA(If J)(J«1.M> 

PRINT  8003 1  (LAC  I ( J) (J=l (H) 

8003  FORMAT ( 1 X > 4013 > 

1  CONTINUE 

C  DEFINE  FREQUENCY  INCREMENTS 

F1=0.5*(FFC2)-FF(1>> 

DO  221  1=1 (NFRQ1 
F2=(FF ( 1+1 )-FF ( I ) )*0<5 
DELF(I)*F1FF2 
221  F1=F2 

DELF ( NFREQ )  =2  • 0*F2 
00  8888  I=1.NFRQ1 

8888  ZFL < I ) =EXP ( - <  FF ( I +1 ) /FF ( I ) ) **4 ) 

C 

DO  921  ITF=1.NFRQ1 

921  DF2(ITF)=(FF(ITF+l)-FF(ITF>>/2. 

C  DEFINE  NONDIHENSIONAL  SHAPE  FUNCTION-FORWARD  FACE-LOCAL  SEA  SOURCE 
C0N1=EXP( 1 • ) 

DO  8849  I=l»95 
FX3=1. 

8849  SNLA(I)=FX3*C0Nl*EXP(-<0.95/<I*0.01>>**4> 

DO  8850  1=95.105 

8850  SNLA( 1 )  =  1  • 

DO  8840  I=104>400 
8840  SNLA( I )=-( 100./I )**5 

C  DEFINE  LIMITING  VALUE  OF  SEA  PEAK  FREQUENCY  AS  FUNCTION  OF  WIND  SPEED 
DX33=(DELX/100. >**(-0.33) 

DO  8855  IU=1>100 
UU=IU*0.514 

F 1=14. 15*UU**( -0.34) *DX33 
8858  V=G/(2.«TW0PI«F1 ) 

V=V*0,01 

D=(DELX/100.  -  V*DELT*0.5)**(-0.33) 

F2=14. 15*UU*«(-0.34)«D 

IF  (ABS(F2-F1).LT. 0.001)  GO  TO  8857 

F1=0.5*(F1+F2) 

GO  TO  8858 

8857  FFLIM< IU)=0.5*(F1+F2) 

8855  FMX( IU) =1 . 74/UU 
18=1 

DO  33  I=l»NFREQ 
II=NFREQ+1-I 
IE=3.04/FF ( 1 1 ) +0 . 5 
DO  34  J=IB. IE 
AFQ7< J)=II 
34  CONTINUE 
IB=IE+1 
33  CONTINUE 

C  DEFINE  NONDIHENSIONAL  SHAPE  FUNCTION  FOR  FORWARD  FACE  OF  SWELL  SOURCE 
DO  8870  1=1 (90 
FX3=(0>90/( 1*0.01 ))**( -3) 

8870  SNLSWC I )=FX3*C0N1*EXP(- (0. 90/( 1*0.01 ) >**4) 

DO  8872  1=91.100 
8872  SNLSWC I )=(95.-I)*0.2 
DO  8875  1=101.400 


A4 


8875  SNLSW(I)*-<100./I)M5 
AINOTWOPI/NDIR 
RADC*TU0PI/360. 

RADDEG=360 . /TUOPI 

C0NEH»GSQ/TW0PI»4 

C0NEE*C0NEH/4. 

EPRH=5.19E-11 

HE0=<FF(NFREQ)+0.5*DELF(NFREQ) )**<-4)*C0NEh/4. 

NF5=NFREQ-5 

ECONST  =6SQ/TU0PIt>4 

NF4=NF5-1 

TPRH=0.1215«T1NC/G*«1. 33333 

EPINRN=2./-PI 

SWOON =9 • B29E-6 

DO  11  KKF=1 *NFREQ 

DO  11  I  =  1*NFREQ. 

FSEA(IiKKF)  =  EXP(-(FF(KKF)/FF(I>)**4) 

11  CONTINUE 

DO  12  I*1*NFRE0 
12  FS£A<1*21)=0. 

DO  651  1=1 i NDIR 
INGLE (I)  =  <I-1)  *  22.5 
ANGL<  I  )  =  (  1-1 XAINC 
COSR ( I ) *COS ( ANGL ( I ) ) 

SINR( I >=SIN( ANGL< I ) ) 

651  CONTINUE 

DO  41  1=1*90 

41  COSSQ(I)  *  COS(d-l)  *  RADC)  **  2 
DO  43  1=91*181 

43  COSSQ(I)  =  0. 

DO  42  1=1*360 
IAAD(l)  =  I 

IF  ( I .GT .180)  IAAD< I )=361-I 

42  CONTINUE 
IAAD1361 )  =  1 

C  ITDIF  IS  ANGLE  DIFF.  BETWEEN  CENTRAL  HAVE  ANGLE  AND  HAVE  ANGLE 
DO  60  ITH=1 *  360 
GAHANGf ITH)=0. 

IF (ITH.GTt80.AND.ITH.LT .280>  GOTO  60 
ANG=(ITH-1)*RADC 

GAHANG<ITH)=C0S(ANG>**4 

GAHANG(ITH)=GANANG(ITH)*8./(3.»PI) 

60  CONTINUE 

C  CONSTANTS  FOR  PHASE  AND  GROUP  VELOCITY 
DO  29  I=1*NFREQ 
Fll ( I )=FF ( I )**1 1 
EFSUM( I )=0. 

C( I )=G/TWOPI/FF<  I ) 

CG( I )=C( I )/2. 

XCR(I)=TINC*CG(I) 

XXD=DELX-XCR(I)/2. 

IF  (XXD.LE . 0. ) . XXD=0. 05 
XD=D£LX/XXD 
RH33 ( I ) =XD**0 . 33 
RH22<I)=XD**0.22 
XCR(I)=XCR<I )/DELX 
29  CONTINUE 

F11<NFREQ+1)=0. 

XCR(NFREQ+1 )=0. 


RM33 ( NFREQ+ 1 ) = 1 . 

RM22  <  MFREQ+ 1 ) ■ 1 . 

DO  70  I«1»NFR£Q 

70  SINF ( I ) »TW0PI*6SQ/ ( TWOPIBFF ( I ) ) **5 
C 

DO  17  IU«1.100 

VST«0.4*31.4*IU 

UST-VST 

19  Z0=C1/UST+C2*UST*UST-C3 
UST1=VST/AL0G( (ZN2-Z0)/Z0> 
IF(ABS(UST1-UST) .LT.0.01 )  GOTO  18 
UST*UST1 

GOTO  19 

18  USTA( IU)*UST1 
C  PRINT  2026.IU.UST1.Z0 
2028  FORMAT  (IX. I5.2F12.6) 

17  CONTINUE 
C 
C 

CALL  CLASSB 
C 

C  BEGIN  NAVE  CALCULATIONS 

C 

C  TEMPORARY  WIND  INPUT 

NUM  =  1000 

IF  (NORD.NE.l)  GO  TO  402 
READ  102.  WINDY. IDIRY.NUM 

102  FORMAT  (313) 

PRINT  103.  WINDY. IDIRY.NUM 

103  FORMAT  (1X.15HTEST  PARAMETERS./. 

.  5X.13HUIND  SPEED  *  >13./. 

.  5X.17HUIND  DIRECTION  «  .13./. 

.  SX.18HNUMBER  OF  HOURS  *  .13) 

DO  401  1=1. N 
DO  401  J=1.M 
WINDA( I  * J)=UINDY 
401  WDIR( I > J)=IDIRY 
402  DO  20  KSS=1 .NPTSP1 
FSWLA(KSS)=1.0 
KFA(KSS)=NFRE0+1 
ESUSQ ( KSS ) =0. 

ISW0(KSS)=NFREQ+1 
IASU(KSS)=0 
DO  8847  NPM=1 .5 
8847  EPMTRX(NPM.XSS)=0, 

FOA(KSS) =0 . 5 

ALPHA(KSS)=0.01 

TH0A(KSS)=0. 

EA(KSS)=1.0E-5 
DO  21  J=1 .NFREO 
EL(KSS» J)=0. 

21  CONTINUE 

DO  20  IA=1 . 16 
DO  20  ITF=1 .20 
TEL(KSS. IA.ITF)  =  0. 

20  CONTINUE 
C 

C  START  CALCULATIONS 
C 


A6 


DO  1000  IHR«1.NUM 
IF(NORD  .EQ.  1)  60  TO  1134 
READ  ( 1 » 1146iEND*9999 )  I  DR 
1146  FORMAT  (110) 

DO  1135  1*1 *N 

READ(01 >1140)  <WINDA(I»J),J»1»M> 

C*t**«*tt 

C  1140  IS  THE  INRUT  FORMAT  FOR  THE  INPUT  WIND  FILE******** 

cm**** 

1140  FORMAT ( 1X>32I3) 

1135  CONTINUE 

DO  1136  1*1 >N 

READ ( 01 > 1140)  <WDIR(IiJ)»J*1fM> 

1136  CONTINUE 
1134  CONTINUE 

9094  DO  1001  KTIMEal »KTIMES 
DO  40  KSS=1 #NPTS 
1SUAR  <  KSS ) =NFREQ+1 
ESUM(KSS)=0. 

USUM(KSS)=0. 

VSUM(KSS)=0. 

AF  07 ( KSS ) *NFREQ+ 1 
FSUM(KSS)=0. 

I-LOCI(KSS) 

J-LOCJ(KSS) 

IUW-UINDA(IfJ) 

USTl-USTA(IUW) 

ALPHA ( KSS )=HBARF (EA(KSS)tUSTl) 

IF ( FO A < KSS ) •  GE>FHX( IUW)  >  GOTO  40 
FOA(KSS)=FMX( IUW) 

KFA(KSS)=AFQ7( IUW) 

40  CONTINUE 
C 

C  DETERMINE  UINDSPEED  AND  DIRECTION  AT  POINT  X-Y 
C 

c***tt***ttttt**t*tttt*tt****tt*t*t***t**t** 

C  PROPAGATION  ROUTINE  FOLLOWED  BY 

C  SOURCE  INTEGRATION 

C  LOOP  STRUCTURE  IS  (FREQ(SPACE(DIRECTION) ) > 

C 

Ctt***t******t**t*******START  OF  FREQUENCY  LOOP  *************** 

C 

C  DEFINE  SOURCE  PARAMETERS 

C 

ISMIN  =  20 

KMIN  =  20 

DO  56  KSS=1  »NPTS. 

IF  (ISWO(KSS).LT. ISMIN)  ISMIN=ISWO(KSS> 

KFRQ=KFA(KSS) 

I=LOCI(KSS) 

J=LOCJ(KSS) 

I UW= W I ND A ( I i J) 

WTH=WDIR( I» J) *RADC 
IANG=WTH/AINC+1 .5 
IFUANG.GE.NDIR)  IANG  =  1 
THO=THOA(KSS) 

IWDIR=TH0A(KSS)/AINC+1.5 
IF(  IUDIR.GT. 16)  IUDIR'l 
IF  (NFRQiLE.NFREQ)  GO  TO  13 


sc*o. 

FQ*1. 

CO  TO  14 

13  IF  <KFRO.LT. KHIN)  KHIN-KFRO 
EHU*XCR<KFRQ) 

IU*I+IUA( IWDIR) 

JU:J+ JUA( IHDIR) 

KSSU'LOCC ( I U  >  JU  > 

IF<KSSU.GT.NPTOT>  GOTO  8343 
EHFU*EHU/2. 

EHF=1.-EHFU 

ALP= ALPHA ( KSS )*EHF+ ALPHA (KSSU)*EHFU 
FO=FOA(KSS)*EHF+FQA<KSSU)*EHFU 
GOTO  8848 
8343  CONTINUE 

IF<FOA(KSS) .LT.FFLIH(IUW))  FOA(KSS)=FFLIH< IUU> 
F0=F0A(KSS)*RH33(KFRQ) 

ALP= ALPHA <  KSS ) KRH22 <  KFRO  > 

8848  SC=2500.«ALP*ALP*ALP/F0**4 

14  SCAA(KSS)=SC 
F000( KSS) =1 . /FO 
TH2=( IANG-1 ) tAINC 
DO  54  IDIR-1 > 16 

IADIF=IABS(WDIR(IrJ>  -  INGLE<IDIR>>  +  1 
I AHT (KSS . IDIR)  =  IAAD< IADIF) 

THDIF=ABS( ANGL( IDIR)-THO) 

ITHDIF=RADDEG*THDIF+1 .01 
ITHDIF=IAAD< ITHDIF) 

IBMT (KSS > IDIR) =ITKDIF 
THDIF=ABS(ANGL(IDIR>-TH2> 

ITHDIF=RADDEGtTHDIF+l ,01 
ITHDIF=IAAB(ITHDIF) 

ICHT(KSS»IDIR)=ITHDIF 
54  CONTINUE 

ISMIN  =  ISHIN  -  5 
KHIN  =  KHIN  -  5 
ITFHIN  =  KHIN 

IF  USHIN.LT. ITFHIN)  ITFHIN  =  ISHIN 
IF  (ITFHIN.LT. 1)  ITFHIN  *  1 
PRINT  8003i ITFHIN 
C 
C 

IF  (ITFHIN.LT. 2)  GO  TO  501 
ITFH1=ITFHIN-1 
DO  500  ITF=1 r ITFM1 
DO  500  KSS=1 »NPTS 

ESUH(KSS)=ESUH(KSS)  +  EL(KSS » ITF >*DELF< ITF ) 

500  CONTINUE 

501  CONTINUE 

DO  50  JJTF=ITFHIN,NFREQ 

I TF= JJTF 

F=FF(ITF) 

F100=F*100. 

KTF=ITF-NF5 
IF (KTF.LT .  1 )  KTF  =  1 
SINFF=SINF(ITF)/AINC 
DO  30  KSS=1  .NPTS 
DO  30  I  A=1 . 16 

E(KSS.IA)  =  TEL(KSS»IA.ITF) 


•  -  -  i 

•  ;■  i 


AS 


30  CONTINUE 
C 

C  PROROGATE  ENERGY 

C 

CALL  PROPR< ITF ) 

C 

C  STRT  OF  SPACE  LOOP. 

C 

DO  57  KSS-1'NPTS 
KFRQ-KFA(KSS) 

EPI=EPHTRX(KTF«KSS> 

SC=SCAA(KSS) 

FOI=FOOO(KSS) 

IFR=F0I*F100 
SCIFR=SC*SNLA< IFR) 

SC2=-2.*SCIFR 

IF  (IFR. LE. 105)  SC2-0. 

IFR=F100/FSULA(KSS) 

SC=SWLSC(KSS) 

SC=SC/(EL(KSS»ITF>  +  0.1) 
SCSU=SCtSNLSU( IFR) 

SUM=0. 

C 

C  START  OF  ANGLE  LOOP 
C 

DO  55  IDIR=liNDIR 
IA=IANT (KSS» IDIR) 

ITHDIF=IBHT(KSSi IDIR) 

JTHDIF=ICHT (KSSi IDIR) 
WNDSRC=SC2*GAMANG< JTHDIF) 
ESP=ENA(KSSrIDIR)+EPI*COSSQ(IA) 
SNL=SCIFRIGANANG< ITHDIF) tSCSUtESP 
EN=ESP+ ( SNL+UNDSRC  >  *T INC 
IF(EN.GE.O. >  GOTO  120 
EN=1.0E-5 

120  EPHON(IDIR)-EN 
SUN=SUN+EN 
55  CONTINUE 
C 

C  END  OF  ANGLE  LOOP 

C 

C  NOTE  SUN=ENERGY/<DELF*AINC) 

C  NQRHAL  SINF  UNITS  ARE  ENERGY/DELF 

C 

SINFKS=ALPHA(NSS)*SINFF 
IF(SUN.LT.SINFKS)  GOTO  68 
RNOR«=SINFKS/SU« 

DO  58  IDIR=1»NDIR 
EPH0N( IDIR) =EPHON( IDIR) fRNORH 
58  CONTINUE 

CALL  DOTPRD(EPHONiCOSRfXSUN) 

CALL  DOTPRD(EPHON»SINR»YSUH) 

USU«(KSS)=USUN(KSS)+XSUH*DELF(ITF> 

VSUN(KSS)=VSUH(KSS)+YSU«*DELF<ITF) 

SUM=SINFKS 

EF ( KSS)=SINFKS 

ESUH(KSS)=ESU«(KSS)  +  SU(UDELF(ITF) 
DO  69  I A=1 t  16 

69  TEL(KSSlIA.ITF)=EPHON(IA) 


GO  TO  57 

68  DO  64  IA*1>  16 

TEL(KSS'IA'ITF)  -  EPHON(IA) 

64  CONTINUE 

AF07(KSS)  -  ITF 
EF<KSS)*SUH 

ESUH(KSS>*ESUH<KSS)+SUM*DELF<ITF) 

57  CONTINUE 

C 

C  END  OF  SPACE  LOOP 

C 

DO  51  KSS-1'NPTS 
51  EL(KSS'ITF)*EF(KSS) 

50  CONTINUE 
C 

C  END  OF  FREQUENCY  LOOP 

C 

DO  66  KSS=1 »NPTS 
HIGHE*HSCALE(KSS) 

I=LOCI(KSS) 

J=LOCJ(KSS) 

IPP(KSS)*NFREQ 
BENG=1 . OE-5 
DO  3007  ITF=1»NFREQ 
EFF=EL<KSS.ITF)*AINC 
IF(EFF.LT.BENG)  GO  TO  3007 
BENG=EFF 
IPP(KSS)=ITF 
3007  CONTINUE 

IU=UINDA(I ' J) 

KFRQ=AF07 ( KSS ) +1 
IMX  =  2 

SNX  *  EL(KSS>2) 

K  *  KFRQ  -  3 
IF  (K.LT.3)  GO  TO  451 
DO  450  ITF  =  3>K 

SEA  =  EL(KSS.ITF)  *  FSEA< ITF'KFRQ) 

SNELL  *  EL < KSS 'ITF)  -  SEA 
IF  < SNELL. LT. SNX)  GO  TO  450 
SNX  *  SNELL 
IMX  *  ITF 

450  CONTINUE 

451  ISN-IMX 

IF  (ISN+3.GE.KFRQ)  GO  TO  22 

ESNL*EL(KSS  > IMX ) +EL(KSS» IHX-1 >4EL(KSS > INX+1 ) 

FS1=EL<KSS'INX)*FF<INX) 

FS2=EL<KSStINX-l)*FF(I»X-l> 

FS3=EL (KSS. IHXF1)*FF< INX+1) 
FSNIA(KSS)*<FS1+FS2+FS3>/ESNL 
ENX*ESHL«0.33*AINC 
ESNSQ<KSS)*EHX*EHX*ENX 

SNLSC ( KSS )*FSNLA< KSS) 1411 tESNSQt  KSS ) 4SNC0N 

EHXA(KSS)*EHX 

GOTO  24 

22  ESNSQCKSS ) *0 • 

FSNLA(KSS)*1 .0 
SWlSC<KSS>*0. 

ISN*NFREQ+1 
24  ISNO(KSS ) = ISW 
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UST1=USTA< IU) 

IF( AF07(KSS)-NFR£Q.LT.0)  GOTO  123 
PARAMETRIC  HIGH  FREOUENCY 
ALP=ALPHA(KSS) 

T=<F0A<KSS>M<-2.33333>+TPRM*<USTl>**<1.33333))**<0. 4285714) 
F0A(KSS)*1./T 
IMDIR=TH0A(KSS>/AINC+1.5 
IF ( IMDIR.GT . 16)  IHDIR=16 
IIU-I4 IUA( IUDIR) 

JU=J+JUA(IUDIR) 

KSSU=L0CC( I IUrJU) 

IF (KSSU.LE.NPTS)  GOTO  7740 
IF(F0A(KSS> .LT.FFLIH(IU))  FOA(KSS) =FFLIM( IU) 

7740  IF(FOA(KSS).LT . FMX( IU) )  FOA(KSS)=FMX< IU) 

C 

C  CALCULATE  ENERGY  INTO  DISCRETE  SPECTRUM  FROM  PARAHETRIC  SOURCE 
C 

EFM=F0A(KSS>**(-5>  *ECONST*EPINRM*ALPHA(KSS) 

SUMH=0. 

DO  8846  NPN=2»5 
NPM=NF4+NPN 

XE=£FM*EXP ( - ( FF ( NPM  > /FOA ( KSS ) ) ** ( -4 . ) ) 

SUMH-SUMH+XEtDELF ( NPM ) 

XED=XE-EPMTRX(NPN»KSS) 

IF1XED.LT.0.)  XED=0, 

8846  EPMTRX(NPN»KSS)=X£D 

THOA(KSS)=UDIR(I.J)*RADC 

KFA(KSS)=NFREQ+1 

FNBAR=F0A<KSS)*UST1/G 

HSCALE  <  KSS  )=EPRM*FMBAR**  <  -3 .333333  >  *UST1«4 
HIGHE=HSCALE(KSS) 

EA<  KSS )  =ESUM(  KSSXAINC+HIGHE 
FSUM(KSS)=FOA<KSS) 

UST1=USTA(IU) 

ALPHA < KSS )=HBARF(EA( KSS) tUSTl ) 

GOTO  66 

123  KFRQ=AF07<KSS) 

KFA(KSS)=KFRQ+1 
IANG=HDIR(I»J>/22.5+1.5 
IF ( IANG.GT .NDIR)  IANG=1 
USUM(KSS)=USUM(KSS)tAINC+HIGHE$COSR( IANG) 

VSUM(KSS)=VSUM(KSS ) *AINC+HIGHE*SINR( I ANG) 

ZMIN=ZFL(KFRQ)*EL(KSSfKFRQ+1 )*AINC 

Z1=EL(KSS.KFRQ)*AINC-ZMIN 

Z2=ALPHA(KSS)*SINF(KFR0)-ZMIN 

Z2-ABS( Z2 ) +0 . 1 

PC=Z1/Z2 

IFIPC.LT. 0. )  PC=0. 

F0A<KSS)=FF(KFRQ+1)-2.*PC*DF2(KFRQ> 

DO  9951  NPN=2.5 
9951  EPMTRX(NPN»KSS)=0. 

FSUM(KSS)=FOA(KSS) 

EA(KSS)=ESUM(KSS) tAINC+HIGHE 
USUM ( KSS )=USUM( KSS 1+0.000001 
X=VSUM(KSS) 

Y=USUM(KSS) 

TH0A(KSS)=ATAN2(X»Y) 


IF (THOA(KSS) .LT.O. )  THOA<KSS)=THOA(KSS>+TUOPI 
UST1=USTA<IU) 

ALPHA <KSS)=HBARF(EA(KSS)  illSTl ) 

HSCAL£(KSS)=HEO* ALPHA (KSS) 

66  CONTINUE 

DO  6660  KSS=1>NPLU 
6660  HLP(KSS)=KFA(KSS)-4 
C 

C  END  OF  TIME  STEP  LOOP 
C 

IF(MSTA.LT.l)  GOTO  4651 
IF(KTIME.NE.l)  GO  TO  1001 
8100  FORMAT < 1H1 f20Xi 10HDAT£-TIME=f 110) 

H0URS=IHR*(TINC*KTIMES)/3600. 

200  FORMAT <5Xf35HTIME  IN  HOURS  SINCE  START  OF  RUN  =  fF6.2f//) 

DO  217  ISTA=1 iHSTA 
PRINT  8100»  IDP 
PRINT  200>  HOURS 
SUM=0. 

I=LOUT ( ISTAi 1 ) 

J=LOUT ( ISTA>2) 

KSS=LOCC(IfJ) 

PRINT  71 r 1 1 J 

71  FORMAT  <3Xf35HSPECIAL  OUTPUT  FOR  GRID  LOCATION 
»2HI=fI3f4H  J  = , 13) 

WRITE<09.8100)  IDP 
DO  8628  ITF=1,20 

URITE(09f8629)  < TEL<KSSi IA» ITF > » IA=1 > 16) 

8628  CONTINUE 

8629  FORMAT ( lXr 16E8<2) 

PRINT  6000 

6000  FORMAT ( 42X . 33HDIMENS IONAL  SPECTRUM  <CM**2/HZ)  ) 

PRINT  6001 

6001  FORMAT  ( 1X> 19HFREQ  ENERGY  DENSITYf25X* 19H  DIRECTIONAL  ENERGY) 
PRINT  6002 

6002  FORMAT ( 1X>4H(HZ) > 2X » 10H(CM<*2/HZ) f30X > 15HDIRECTI0N  BANDS) 

PRINT  199 

199  FORMAT (19Xf,0.0'f4Xf,22*5*f3Xf,45»0,f3Xf,67.5,f3Xf,90.0,f 

*2Xf* 112.5* f2Xf' 135.0' f2Xf‘ 157. 5* f2Xf ' 180.0* f2Xf *202.5' f2Xf *225.0 
*,f2Xf’247.5’f2Xf,270.0,f2Xf*292.5,f2Xf,315.0,f2Xf,337.5,f//) 

C 

C  INTEGRATE  TO  GET  ONE-DIMENSIONAL  SPECTRUM 
C 

WRITE(08f8100)  IDP 
FSS=0.0 

DO  317  ITF=1fNFREQ 

EFS=EL(KSSfITF)*AINC 

SUM=SUM+EFS*DELF(ITF) 

FSS=FSS+EFS*DELF( ITF)*FF( ITF ) 

PRINT  73 f  FF(ITF)fEFSf< TEL<KSSf IAf ITF) f IA=1 r 16) 

73  FORMAT  ( 1XfF5. 3f 1XfF7 . Of IXf 16F7 . 0 > 

HRITE(08f76)  ITFfEFSfDELF<ITF)fSUM 
76  F0RMAT(1XfI3.3E13.6) 

317  CONTINUE 

PP0SP=1.0/FF(IPP(KSS)) 

C  PPOSP=  PEAK  OF  SPECTRUM  IN  DISCREET  REGION 

IF(FSS.LE.O.O)  FSS=1 . OE-5 
TSS-SUM/FSS 
C 


C  PRINT  ONE-  AND  TWO-DIMENSIONAL  SPECTRUM 
C 

PRINT  6004, SUM, HSC ALE (KSS) 

6004  FORMAT  <  5X ,  1 2HT0T  AL  ENERGY, Ell. 4, 3H  +  *  El 1 . 4, 13H( PARAMETRIC)  ) 
SUM=SUM+HSCALE(KSS) 

USTl=USTA(WINDA(I,J)) 

EBAR=GSQ*SUM/UST1**4 
HBAR=SQRT (EBAR) 

HT=4.01*SQRT(SUM> 

HTT=.01*HT 
PRINT  6005, HTT 

6005  FORMAT <5X,34HSIGNIFICANT  WAVE  HEIGHT  (METERS)  =,9X»F5.1) 
HT2=HT*0,0328 

PRINT  315,  HT2 

315  FORMAT  <5X»32HSIGNIFICANT  WAVE  HEIGHT  (FEET)  =,11X,F5.1) 
TPER=1 .0/FQA<  KSS) 

PRINT  6006, TPER 

6006  FORMAT (5X,34HPERI0D  OF  THE  SEA  PEAK  (SECONDS)  =>9X,F5.1) 

PRINT  6020, TSS 

6020  FORMA T(5X»39HAVERAGE  PERIOD  OF  ALL  WAVES  (SECONDS)  =,4X,F5.1) 

PRINT  6017»PP0SP 

6017  FORMAT (5X,35HPEAK  PERIOD  OF  SPECTRUM  (SECONDS)  =,8X,F5.1) 

ITH0S=TH2*RADDEG 
PRINT  6007, ITHOS 

6007  FORMAT (5X»37HME AN  DIRECTION  OF  THE  SEA(DEGREES)  =  ,6X,IS) 
ITHOP=THOA(KSS)*RADDEG 

PRINT  6008, ITHOP 

6008  FORMAT (5X,38HMEAN  DIRECTION  OF  ALL  WAVES  (DEGREES)1 > 5X> 15) 
PRINT  6009,WINDA( I , J) 

6009  FORMAT (5X,34HWIND  SPEED  AT  GRID  POINT(KNOTS)  =,9X,I5> 

PRINT  6010,WDIR( I , J) 

6010  FORMAT (5X,43HDIRECTI0N  OF  WIND  AT  GRID  POINT (DEGREES)=  ,15) 
WRITE(08,365)  HT ,HT2 

WRITE(08,316)HBAR»EBAR 

316  FORMAT  (1X,2E12.4) 

365  FORMAT ( 1X,2F10.2) 

217  CONTINUE 

C 

C  PRINT  INFORMATION  FOR  ENTIRE  GRID 
C 

4651  DO  8050  1=1, N 
DO  8050  J=1,M 
H2(r,J)=0. 

8050  CONTINUE 

DO  444  KSS=1,NPTS 
I=LOCI (KSS ) 

J=LOCJ(KSS) 

444  H2(I,J)=0.0401*SQRT(EA(KSS>> 

6666  FORMAT  (1H1) 

PRINT  8100,  1DP 
PRINT  6011 

6011  FORMAT  (5X,47HSPATIAL  GRID-SIGNIFICANT  WAVE  HEIGHT (METERS)  ) 
DO  311  1  =  1,  N 

PRINT  312,  <H2(I,J),J=1,M> 

312  FORMAT ( IX, 15F6. 1 ) 

311  CONTINUE 

DO  313  KSS=1,NPTS 
I=LOCI (KSS) 

J=LOCJ(KSS) 


313  H2( I . J)=THOA(KSS) (RADDEGtO. 1 

PRINT  4012 

6012  FORMAT (5X.51HSPATIAL  GRID-MEAN  DIRECTION  OF  SPECTRUM(DEGREES/10>  ) 
DO  445  1=1 tN 

445  PRINT  312*  (H2(I. J) . J=1.M) 

PRINT  9100 

9100  FORMAT (5X*33HSPATIAL  GRID-PEAK  PERIOD(SECONDS) ) 

DO  314  KSS=1*NPTS 
I=LOCI ( KSS) 

J=LOCJ(KSS) 

314  H2( I » J)  =  l .0/FF( IPP(KSS)  > 

DO  446  1=1 *N 

446  PRINT  312.  (H2< I  * J>  * J=1 *M) 

1001  CONTINUE 

1000  CONTINUE 
9999  STOP 
END 

SUBROUTINE  CLASSB 

COMMON  /DD/CG(20) .DELXiTINC 

COMMON  /ST/XMA1 (16*20)*XMA2(16>20) ,PU1 (20)  *PU2(  20) 

COMMON  /CB/LA(31.31)fL0CC<36.36) 

COMMON  /CC/  N.M.NPLU.NPNLy.NPTS.NFREQ.NDlR.NPTSPl 
COMMON  /EE/LOCI (961 ) *L0CJ<961 ) 

COMMON  /RP/KSSULU( 729. 16) .EV(5. 16) .LXLI ( 16.5) .LXLJ( 16.5). 

1 1 REF  A  < 16) 

COMMON  /KU/  KSSUA (232. 16) 

COMMON  /IJ/  IUA ( 16) . JUA( 16 ) 

COMMON  /B/  VMU(5. 16.20 ) .PMU( 5.20) 

KSS=0 

NN=N-2 

MM=M-2 

C  DETERMINE  FULLY  LAX-UENDROFF  REGION 
DO  1  1=3. NN 
DO  1  J=3 . MM 
LOCC ( I . J) =0 

IF(LAd.J).LE.O)  GOTO  1 
IF(LA(I+l.J),LE.O)  GOTO  1 
IF(LA(I+l,J+l).LE.O)  GOTO  1 
IF ( LA ( I . J+l ) ,LE . 0)  GOTO  1 
IF(LA(I-1.J+1).LE.0>  GOTO  1 
IF ( LA( I -1 > J )  ,LE . 0 )  GOTO  1 
IF(LAd-l.J-l).LE.O)  GOTO  1 
IF  ( LA ( I .  J- 1 )  .LE .  0 )  GOTO  1 
IF <LA( 1+1 >J— 1) .LE .0)  GOTO  1 
IF ( LA ( I +2. J ) . LE . 0 )  GOTO  1 
IF(LA(I+2.JF2).LE.O)  GOTO  1 
IF(LA(I,J+2).LE.O)  GOTO  1 
IF(LA(I-2.J+2).LE.0>  GOTO  1 
IF(LA( 1-2* J) .LE.O)  GOTO  1 
IF (LA( 1-2. J-2) . LE.O)  GOTO  1 
I F < L A ( I . J-2) .LE.O)  GOTO  1 
IF ( LA( 1+2. J-2). LE.O)  GOTO  1 
KSS=KSS+1 
LOCC(I.J)=KSS 
LOCI (KSS)=I 
LOCJ(KSS)=J 
1  CONTINUE 
NF'LM=KSS 

C  NPLU  IS  THE  NUMBER  OF  POINTS  IN  THE  LAX-NENDROFF  REGION 


4 


A14 


1  I 


NPNLU=0 
N1=N-1 
M1=M-1 
00  2  1=2. HI 
DO  2  J=2»H1 

IF (LOCC ( I , J) . GE . 1 )  GOTO  2 

IF (LA( I , J ) .L£ • 0)  GOTO  2 

NPNLM=NPNLU+1 

KSS=KSS+1 

LOCC ( I i J ) =KSS 

LOCI ( KSS )  =  1 

LQCJ(KSS)=J 

2  CONTINUE 

NPNLU  IS  NUMBER  OF  BOUNDARY  REGION  POINTS 
KSS  NOW  EQUALS  TOTAL  POINTS 
DO  500  1  =  1.  N 

PRINT  501.  (LOCC ( I . J) . J=1 ,H) 

501  FORMAT ( IX. 3114) 

500  CONTINUE 
NPTS=KSS 
NPTSP1=NPTS+1 

PRINT  503.  NPTS.NPLU.NPNLU 
503  F0RMAT<5X. 31 10) 

DO  3  1=1. N 
DO  3  J=l»h 

IF ( LA ( I » J ) . EQ . 0 )  LOCC ( I » J)=NPTSP1 

3  CONTINUE 

CALCULATE  AND  STORE  MULTIPLIERS  FOR  BOUNDARY  REGION 
DO  40  ITF=1 .NFREQ 
DIST=CG( ITF )*TINC 
D22=0.9229*DIST 
D45=0.707*DIST 
DO  4  IA=1 . 13,4 
XMA2( I  A, ITF)=DIST/DELX 
XMA1(IA,ITF)=1.0-XMA2(IA»ITF) 

4  CONTINUE 

DO  5  I A=3 .15.4 

XMA2( IA, ITF)=D45/DELX 

XMA1  ( IA» ITF)  =  1 .0  *.  42(IArITF> 

5  CONTINUE 

DO  6  I A=2 ,16.2 

XMA2(IA.ITF)=D22/DELX 

XMA1(IA.ITF)=1.0-XMA2(IA»ITF) 

6  CONTINUE 

DO  7  I A= 1 , 16 
ENU=XMA2< IA, ITF) 

EMUSQ=EMU*EMU 
EMUP=1 . +EMU 
EMU2=EMU/2. 

EMUM=  1 . -EMU 
A=l.+3.*(l.-EMUSQ)/4. 

B  =  < 1 . -EMUSQ)/4 . 

UMU(1,IA,ITF)=1.-EMUSQ*A 

VMU(2,IA,ITF)=EMU2*A*EMUP-EMU2*B*EMUM 

VMU(3,IA,ITF)=EMU2»B*EMUP 

VMU<  4, IA, ITF) =EMU2*( BIEMUP-A4EMUM) 

VMU(5,IA,ITF)=B*EMU2*EMUM 

7  CONTINUE 

PU2(ITF)=0.3827*DIST/DELX 


PU1 ( ITF)=1 •-PU2( ITF) 

EMU=PU2( ITF) 

EMUSQ=EMU*EHU 

EHUP=1.0+EHU 

EMU2=EMU/2. 

EMUM=1.-EMU 

A=l.+3.*(l.-ENUSQ)/4. 

B=(l.-EMUSQ)/4. 

PHU( 1 1 ITF  > =1 . -EMUSQt A 

PMU(2fITF)=EMU2*A*EMUP-EMU2*B*EMUM 

PMU(3fITF)=EMU2*B*EHUP 

PMU(4f ITF )=EHU2><  B*£MUP-A<EHUH) 

PMU(5f ITF) =B*£MU2*ENUM 
40  CONTINUE 

DO  7716  K=1fNPNLW 
KSS=K+NPLU 
I-LOCI (KSS) 

J-LOCJ(KSS) 

DO  7716  IA=1fNDIR 
1 1  =  I A 

I U= I + IUA (II) 

JU=J+JUA(II) 

KSSUA(Kf IA)=LOCC( IUi JU) 

7716  CONTINUE 
NPTSP1=NPTS+1 

DO  7717  KSS-1 rNPLU 
I=LOCI (KSS) 

J=LOCJ(KSS) 

DO  7717  IA=1 r 16 
I U= I + I UA ( IA) 

JU=J+JUA(IA) 

KSSULU ( KSS . I A ) =LOCC  < IU » JU) 

7717  CONTINUE 

DO  7718  K=2f5 

KK=K 

KH=1 

IF(M0D(KKf2>.EQ.1)  K«=2 
DO  7718  IA=1  f 16 
LXLI(IAfKK)=IUA(IA)*KM 
LXLJ(IAfKK)*JUA(IA)*KN 

7718  CONTINUE 
RETURN 
END 

SUBROUTINE  PROPR  <ITF> 

COMMON  /A/  E(961f16)fEN(961f16) 

COMMON  /ST/XMA1  <  1 6 f 20 )  fXMA2(  16f20)fPU1(20)fPU2(20) 

COMMON  /CB/LA(31f31)fL0CC(36f36) 

COMMON  /CC/  NfMfNPLNfNPNLUfNPTSfNFREQfNDIRfNPTSPI 
COMMON  /RP/KSSULU< 729f 16)fEV(5f16)fLXLI(56f5)fLXLJ(16f5) 
. fIREFA(16) 

COMMON  /»/  VMU(Sf16f20)fPMU(5f20) 

COMMON  /EE /LOCI (961 ) fLOC J (961 > 

COMMON  /KU/KSSUA(232f 16) 

COMMON  /NP/KFA(961 ) fMLP(961  ) 

NTOT=NPTS 

DO  1  KSS=1fNPLU 

IF  (ITF.LT.MLP(KSS))  GO  TO  2 

DO  4  I A= 1 f 16 

KSSU=KSSULU(KSSfIA) 


Ala 


ENIKSSf IA )>XMA1  ( I Af ITF ) IE ( KSSf I A ) +XHA2 ( I A f ITF  >*E ( KSSUf  IA) 

4  CONTINUE 
GO  TO  1 

2  I-LOCI (KSS ) 

J*LOCJ(KSS) 

DO  3  IA=lrl6 
EV(1»IA)=E(KSS.IA> 

DO  44  K=2»5 
LI=I+LXLI < IAiK) 

LJ=J+LXLJ(IA»K> 

KSL=LOCC(LI»LJ) 

EV(Kr IA)=E(KSLfIA) 

44  CONTINUE 

3  CONTINUE 

DO  6  I A=1 , 16 

ENIKSSf IA) -EV( 1 1 IA) tVMU< 1 f IA» ITF)+ 

.  EV(  2 f  I AXUHU(  2f  1  At  ITF )  + 

.  EV(3f  IAXVHU(3f  IAr  ITF )  + 

.  EU<  4f  I A  XVHU(4  f  I  Af ITF)  + 

.  EV(5fIAXVHU(5fIAfITF) 

6  CONTINUE 

I  CONTINUE 

DO  190  K=1 fNPNLN 
KSS=K+NPLU 
DO  141  IA=1 » 16 
KSSU=KSSUA(K> I  A) 

EN(  KSSi  IA)=XNA1  ( IAf  ITF  XE(KSS  f  I A ) +XMA2  ( I  Af  ITF  XE  (  KSSUf  IA ) 
141  CONTINUE 
190  CONTINUE 

C  2ND  HALF  OF  PROPAGATION  IN  LAX-UENDROFF  REGION 
DO  11  KSS=1 fNPLU 
IF  < ITF.LE.NLP(KSS) )  GO  TO  12 
DO  5  IA=2f16f2 
IIA=IREFA< IA) 

KSSUsKSSULN(KSSfIIA) 

E ( KSS . I  A) =PU1 ( ITF ) *EN(KSS f I  A) +PU2( ITF ) »EN( KSSUf I A) 

5  CONTINUE 
GO  TO  11 

12  DO  13  IA=2f 16f  2 
IIA=IREFA< IA) 

EV(1fIA)=EN(KSSfIA) 

DO  14  K=2f5 
LI=I+LXLI<IIAfK) 

LJ=J+LXLJ(IIAfK> 

KSL=LOCC<LIfLJ> 

EV(KfIA)=EN(KSLfIA) 

14  CONTINUE 

13  CONTINUE 

DO  16  IA=2f 16 f 2 
E(KSSfIA)=EV<1fIA)*PNU(1fITF) 

.  +EV(2fIAXPMU(2fITF> 

.  +EV<  3f  IAXPHU(3f  ITF) 

.  +EV<4fIAXPMU<4fITF> 

.  +EV(5  f I AXPHU  (5  f ITF ) 

16  CONTINUE 

II  CONTINUE 

C  2ND  HALF  OF  PROPAGATION  IN  BOUNDARY  REGION 
DO  9  K=1fNPNLW 
KSS=KFNPLW 


00  41  IA«2*16*2 
IIA«IREFA(IA) 

KSSU’KSSUA ( K  *  1 1 A ) 

E (KSS> I A) *PU1 ( ITF )*£N( KSSt I A ) +PU2( ITF )*EN(KSSU> I A) 

CONTINUE 

CONTINUE 

DO  7  KSS=1 *NTOT 

DO  7  IA»2> 16.2 

EN<KSS*IA)*E<KSSiIA> 

CONTINUE 
DO  8  IA=1*16 
DO  8  KSS*1 tNTOT 

IF  <£N(KSS* IA) . LT. 0. )  EN(KSS* I A ) =0 • 

CONTINUE 

RETURN 

END 

FUNCTION  HBARF(EfUST) 

UST4=UST**4 

EBAR -960400 . tE/UST  4 

HBARF=0.044*EBAR**(-0.2> 

IF  (HBARF.LT. 0.008)  HBARF=0.008 

RETURN 

END 

SUBROUTINE  DOTPRD  <X,Y,Z) 

DIMENSION  X ( 16) »  Y < 16 ) 

Z=0. 

DO  1  1=1*16 

z=z+x<n*r(i) 

1  CONTINUE 
RETURN 
END 


APPENDIX  B:  MODIFICATIONS  FOR  A  SPHERICAL  ORTHOGONAL  GRID 


1.  The  flat  earth  grid  will  give  good  results  for  small  bodies 
of  water;  but  if  a  large  body  of  water  is  being  considered,  a  modifi¬ 
cation  to  represent  the  grid  as  a  section  of  a  spherical  earth  is 
needed.  The  spherical  orthogonal  grid  requires  changes  to  be  made  in 
the  dimension  statements  and  in  the  main  and  subroutine  portions  of  the 
program.  A  listing  of  the  changed  subroutines  appears  at  the  end  of 
this  appendix.  This  change  will  cause  the  core  storage  to  increase 
significantly.  The  spherical  orthogonal  grid  requires  that  the  grid 
blocks  change  from  the  nearly  square  blocks  at  the  center  of  the  grid 
to  rectangular  blocks  at  the  edge  of  the  grid.  Changes  involve  an 
additional  array  for  DELY,  the  distance  between  grid  points  in  the  y- 
direction.  The  y-distances  will  differ  depending  on  the  J  counter  (J 
is  the  grid  column  counter)  in  the  program.  DELX  will  remain  the  same 
at  all  J  values.  DELY  (dimensioned  by  the  J  counter)  can  be  added  to 
the  common  /DD/  block  in  the  dimension  section.  In  order  to  define 
DELY  two  new  variables,  JCEN  and  SDEG  (J-center  of  the  grid  and  grid 
spacing  in  degrees),  must  be  read  into  the  program.  These  values  de¬ 
pend  on  the  specific  grid.  The  following  loop  will  set  up  the  DELY 
distances  for  each  J  counter. 

DO  659  J  =  1,M 

DIF  =  ABS(J-JCEN)*SDEG*RADC 

DELY(J)  =  DELX*COS (DIF) 

659  CONTINUE 

'  >e  J  grid  values  go  over  the  entire  grid  where  M  is  the  maximum  J. 

RADC  is  the  conversion  from  degrees  to  radians.  This  loop  needs  to  be 
inserted  early  in  the  main  program  before  subroutine  CLASSB  is  called. 

2.  Other  dimension  changes  in  the  common  block  statements  include 
the  variables,  XMA1,  XMA2,  PU1,  PU2 ,  VMU,  and  PMU.  All  these  variables 
need  to  have  a  dimension  added  so  each  of  them  will  allow  different 
values  at  each  different  J  counter.  These  variables  will  now  be  three- 
and  four-dimensional  variables  and  core  storage  will  increase  by  a 
factor  of  the  maximum  value  of  J  times  the  initial  storage  value  for 


these  variables.  These  variables  are  used  in  the  two  large  sub¬ 
routines,  CLASSB  and  PROPR. 

3.  The  subroutine  CLASSB  sets  up  the  constants  that  are  used  in 
the  propagation  routine.  These  constants  need  to  be  changed  to  fit  the 
varying  block  size  needed  by  the  spherical  orthogonal  changes.  The  DO 
loops  referring  to  statement  numbers  4,  5,  and  6  calculate  the  distance 
values  of  the  block  in  the  various  directions.  These  loops  correspond 
to  IA  =  1(0°),  IA  =  2(22.5°),  IA  =  3(45°),  and  so  on  up  to  IA  =  16 
(337.5°). 

4.  The  DO  4  loop  calculates  the  distances  for  0°,  90°,  180®, 
and  270°.  Change  DELX  to  DELY(J).  If  an  IF  statement 

IF  (IA.EQ.1.0R.9)  DELY(J)  =  DELX 

is  added  after  the  DO  4,  DELX  will  be  used  instead  of  DELY  for  the  0° 
and  180°  distance  increment.  A  second  DO  4  index 
DO  4  J  =  1 ,M 

needs  to  be  added  to  the  loop  to  calculate  all  the  different  DELY(J) 
values.  Obviously,  XMA2  and  XMA1  need  to  have  a  J  dimension  added. 

5.  The  DO  loop  corresponding  to  the  5  statement  number  in 
CLASSB  calculates  the  propagation  constants  in  the  IA  =  3(45°), 

IA  =  7(135°),  IA  =  11(225°),  and  IA  =  15(315°).  DELX  needs  to  be 
changed  to  DELY(J).  Again,  XMA2  and  XMA1  need  to  be  dimensioned  by 
J  and  a 

DO  5  J  =  1,M 

needs  to  be  added  at  the  beginning  of  the  loop  to  index  all  the  J 
values . 

6.  The  DO  6  loop  calculates  the  22.5°,  67.5°,  112.5°,  157.5°, 
202.5°,  247.5°,  292.5°,  and  337.5°  propagation  constants.  DELX  needs 
to  be  changed  to  DELY(J).  An  IF  statement  needs  to  be  added  in  the 
loop: 

IF  (IA.EQ.2.0R.8.0R. 10.0R. 16)  DELY(J)  =  DELX 
As  before,  XMA2  and  XMA1  need  to  be  dimensioned  by  J  and 
DO  6  J  =  1  ,M 

needs  to  be  added  at  the  beginning  of  the  loop  to  index  all  the  J 


7.  The  DO  7  loop  calculates  values  for  the  VMU  variable.  This 
loop  needs  to  have  a 

DO  7  J  =  1,11 

added  to  index  the  various  J  values .  All  the  VMU  variables  need  to  be 
indexed  and  dimensioned  by  J. 

8.  A  loop  needs  to  be  added  for  the  calculation  of  the  PU2  and 
PU1  variables.  These  variables  involve  the  second  component  (the  other 
component  was  calculated  previously)  of  the  propagation  constant  for 
the  22.5°  and  67.5°  angle  rays.  A  suggestion  for  this  loop  is 

DO  AO  J=1,M 
DO  40  IA=2, 16,2 

IF  (IA.EQ.4.0R.6.0R. 12.0R. 14)  DELY(J)=DELX 
PU2  (IA, J,ITF)=0.3827*DIST/DELY(J) 

PU1  (IA, J, ITF)=1-PU2  (IA, J, ITF) 

EMU=PU2  (IA, J, ITF) 

EMUSQ=EMU*EMU 
EMUP=1 .0+EMU 
EMU2=EMU/2 . 

EMUM=1 . -EMU 
A=l.+3.*(l.-EMUSQ)/4. 

B=(l.-EMUSQ)/4. 

PMU ( 1 , I A , J , ITF ) = 1 . -EMUSQ*A 

PMU(2 , IA, J, ITF)=EMU2*A*EMUP-EMU2*B*EMUM 

PMU ( 3 , I A , J , ITF ) =-EMU2*B*EMUP 

PMU(4 , IA , J , ITF)=EMU2*(B*EMUP-A*EMUM) 

PMU (5 , I A , J , ITF) =B*EMU2*EMUM 

40  CONTINUE 

The  dimensioning  changes  of  the  PMU  variables  are  also  shown  in  the 
above  loop. 

9.  The  changes  in  the  propagation  subroutine  PROPR  are  rela¬ 
tively  easy.  The  variables  just  discussed  in  CLASSB  need  to  be 
redimensioned  in  PROPR.  The  J  dimension  for  the  grid  point  (KSS 
point)  needs  to  be  known  to  utilize  the  correct  propagation  constants. 
JKSS=LOCJ(KSS)  gives  the  needed  J  dimension  and  needs  to  be  inserted 
in  the  various  loops  to  assure  that  the  proper  propagation  constants 
are  being  used. 

10.  A  listing  of  the  two  subroutines  that  have  been  modified  for 
the  spherical  orthogonal  grid  follows.  Note  that  the  loop  discussed  in 
paragraph  1  must  be  included  in  the  main  program  and  the  dimension 


statements  in  the  main  program  must  correspond  to  the  dimension 
statements  in  the  subroutines.  The  main  program  must  also  receive 

SDEG  and  JCEN  as  described  in  paragraph  1  as  input. 

SUBROUTINE  CLASSB 

COMMON  /DD/CG( 20 ) >  DELX>  TINC  >  DELY ( 31 > 

COMMON  /ST/XMA1 ( 16>31 >20) >XMA2( 16>31 >20) >PU1 < 16>  31 >  20)  > 

1  PU2( 16  >31 > 20) 

COMMON  /CB/LA<31 >31 ) >L0CC(36>36) 

COMMON  /CC/  N>M>NPLU>NPNLU>NPTS>NFREQ>NDIR>NPTSP1 
COMMON  /EE/LOCI <961 ) >L0C J(961 > 

COMMON  /RP/KSSULU <  729> 16) >  EV( 5  > 16) fLXLI ( 16 > 5 ) >LXLJ( 16>3)  > 
1IREFAU6) 

COMMON  /KU/  KSSUA t 232> 16) 

COMMON  /IJ/  I U A ( 16) > JUA(16> 

COMMON  /B/  UHU(S> 16>3i >20) >PMU(S> 16>31 >20) 

KSS=0 

NN=N-2 

MM=M-2 

C  DETERMINE  FULLY  LAX-UENDROFF  REGION 
DO  1  1-2, NH 
DO  1  J=3>MM 
L0CC(I> J)-0 

IF ( LA ( I > J ) . LE . 0 )  GOTO  1 
IF ( LAI  1  +  1 >  J ) .LE .0)  GOTO  1 
IF<LA(I+l>J+l).LE.O)  GOTO  1 
IF(LA( I > J+l ) .LE.O)  OOTO  1 
IF(LA(I-1>J+1).LE.0)  GOTO  1 
IF(LA( 1-1 > J) .LE.O)  GOTO  1 
IF(LA(I-l>J-l).LE.O)  GOTO  1 
IF ( LA ( I > J-l > .LE .0)  GOTO  1 
IF(LA(I+l>J-l).LE.O)  GOTO  1 
IF<LA<I+2>J).LE.O)  GOTO  1 
IF<LA<I+2>J+2). LE.O)  GOTO  1 
I F ( LA < I » J  +  2 ) .LE.O)  GOTO  1 
IF <LA( I-2> J+2) .LE.O)  GOTO  1 
IF ( LA ( 1-2 > J ) . LE .0 )  GOTO  1 
IF(LA( I-2> J-2) .LE. 0)  GOTO  1 
IF(LA( I > J-2) .LE.O)  GOTO  1 
IF(LA( I+2> J-2) .LE. 0)  GOTO  1 
KSS=KSS+1 
L0CC( I > J)-KSS 
LOCI (KSS)=I 
L0CJ(KSS)=J 

1  CONTINUE 
NPLU=KSS 

C  NPLU  IS  THE  NUMBER  OF  POINTS  IN  THE  LAX-UENDROFF  REGION 
NPNLU=0 
N1=N-1 
M1=M-1 
DO  2  I=2>N1 
DO  2  J=2>M1 

IF (L0CC( I > J) .GE. 1 )  GOTO  2 

IF(LA(I > J ) .LE.O )  GOTO  2 

NPNLU=NPNLU+1 

KSS=KSS+1 

LOCC ( I >  J) =KSS 

LOCI (KSS) =1 

L0CJ(KSS)=J 

2  CONTINUE 

C  NPNLU  IS  NUMBER  OF  BOUNDARY  REGION  POINTS 
C  KSS  NOU  EQUALS  TOTAL  POINTS 


DO  500  I  =  1»N 

PRINT  501  f  ( LOCC ( I »J)iJ=lrM) 

501  FORMAT (lXr 3114) 

500  CONTINUE 
NPTS=KSS 
NPTSP1=NPTS+ 1 

PRINT  503 1  NPTS»NPLH»NPNLH 
503  F0RMAT(5Xi 3110) 

DO  3  1=1 rN 
DO  3  J=liM 

IF(LAdfJ).LE.O)  LOCC(  I r J)=NPTSP1 
3  CONTINUE 

CALCULATE  AND  STORE  MULTIPLIERS  FOR  BOUNDARY  REGION 
DO  40  ITF=1 1 NFREQ 
DIST=CG( ITF) tTINC 
D22=0.9239*DIST 
D45=0.707*DIST 
DO  4  J=1.H 
DO  4  1A=1i13>4 

IF  (IA.EQ.1.0R.9)  DELY(J)=DELX 
XNA2( IAi Ji ITF )=DIST/DELY( J) 

XMA1 < IAi Ji ITF )  =  1 .0-XMA21 1  At J» ITF ) 

4  CONTINUE 

DO  5  J=liM 
DO  5  I A=3 i 15>4 
XHA2( IAi J» ITF) =D45/D£LY ( J) 

XNA1(IA*J.ITF)=1 .0-XNA2< I A» Ji ITF) 

5  CONTINUE 

DO  6  J=1 »N 
DO  6  1A=2.16'2 

IFUA.EQ.2.0R.8.0R.10.0R.16)  DELY(J)=DELX 
XMA2( IAi Ji ITF )=D22/DELY( J) 

XMA1 < IAi  Ji ITF )  =  1 .0-XMA2( IA» Ji ITF) 

6  CONTINUE 

DO  7  J=1.N 
DO  7  I A = 1  *  1 6 
EMU=XNA2( IAi J» ITF) 

EMUSQ=EMUtEHU 
EMUP=1 . +ENU 
E«U2=E«U/2. 

ENUH=1 .-EMU 
A=l.+3.*(l.-EMUSQ)/4. 

B=(l.-EMUSQ)/4. 

VMUUiIAi  JiITF)  =  l.-ENUSQ*A 
VMU(2rIA»JrITF)=EMU2*A*EMUP-EMU2*B*EMUM 
VMU<3iIAi JiITF)=-ENU2*B*EMUP 
VMU<4,IAf JfITF)=EMU2*(B*EMUP-A*EMUM) 
VMU(5iIA.JiITF)=B*EMU2*EMUH 

7  CONTINUE 

DO  40  J=1 »M 

DO  40  IA=2i 16 12 

IFUA.EQ.4.0R.6.0R.12.0R.14)  DELY(J)=DELX 
PU2( IAi Ji ITF  >=0. 38274DIST/DELY ( J) 

PUKIAf  JfITF)  =  l.-PU2(IAf  JfITF) 

EMU=PU2< IAi Ji ITF) 

EMUSQ=EHU*EMU 
EMUP=1 .O+EMU 
EMU2=EMU/2. 

EMUM=1 . -EMU 


A*l.+3.*<l.-EMUSQ)/4. 

B»(l.-EMUS0>/4. 

PhU  < 1 fIAfJfITF)«1. -EMUSQSA 

PHU<2»  IAt  J»ITF )  =EMU2*A*EHUP-EMU2*B*EMUM 

PMU(3fIAfJ,ITF)*-EMU2*B*EMUP 

PMU(4f IAf Jf ITF ) =EMU2t(B*EMUP-A*EMUM> 

PNU  <  5 » I A .  J .  I TF )  *B*EMU2*EHUM 
40  CONTINUE 

DO  7716  KMfNPNLU 
KSS=K+NPLH 
I=LOCI (KSS) 

J=LOCJ(KSS) 

DO  7716  IA=1fNDIR 
II  =  IA 

IU=I  +  IUA( II ) 

JU  =  J+ JUA( II ) 

KSSUA<KfIA)*LOCC(IUfJU) 

7716  CONTINUE 
NPTSP1=NPTS+1 

DO  7717  KSS=1fNPLU 
I=LOCI (KSS) 

J-LOCJ(KSS) 

DO  7717  IA=1 f 16 
IU=I+IUA(IA) 

JU=J+JUA( IA) 

KSSULU(KSSfIA)=LOCC(IUfJU) 

7717  CONTINUE 

DO  7718  K=2f5 

KK=K 

KM=1 

IF(«0D(KK.2).EQ.l)  KM=2 
DO  7718  IA=1f16 
LXLI<IAfKK)=IUA(IA)»KH 
LXLJ(IAfKK)=JUA<IA)*KH 

7718  CONTINUE 
RETURN 
END 

SUBROUTINE  PROPR  (ITF) 

COMMON  /A/  E(961f16)fEN<961f16> 

COMMON  /ST/XMA1 ( 1 6 f 3 1 f 20 > fXHA2< 16f31f20>fPU1(16f31f20)f 
1  PU2( 16f31 f 20) 

COMMON  /CB/LA( 31 f31 ) fLOCC ( 36  f  36  > 

COMMON  /CC/  NfMfNPLWfNPNLUfNPTSfNFREQfNDIRfNPTSPI 
COMMON  /RP/KSSULU(729f16)fEM(5f16)fLXLI(16f5)fLXLJ(16f5) 

. fIREFA( 16) 

COMMON  /B/  UMU(5f16f31f20) fPMU(5f 16 f 31 f 20) 

COMMON  /EE/LOCI <961 ) fLOCJI 961 ) 

COMMON  /KU/KSSUA(232i 16) 

COMMON  /NP/KFA(962)fMLP(961) 

NTOT=NPTS 
DO  1  KSS=1fNPLU 
JKSS=LOCJ(KSS) 

IF  ( ITF. IT .MLP(KSS) )  GO  TO  2 
DO  4  I A-l f  16 
KSSU=KSSULH(KSSfIA) 

EN(KSSfIA)=XMA1(IAfJKSSfITF)*E(KSSfIA)+XMA2(IAfJKSSfITF)*E(KSSU 
l  F  IA) 

4  CONTINUE 
GO  TO  1 


2 


I’LOCKKSS) 

J-LOCJ(KSS) 

DO  3  IA=1f16 
EV< 1 f IA)*E(KSSf IA) 

DO  44  K=2»5 
LI*I+LXLI(IAfK) 

L J=J+LXLJ( IAfK) 

KSL=LOCC(LI >LJ) 

EV<KfIA)*E(KSLfIA> 

44  CONTINUE 
3  CONTINUE 

DO  6  IA-1 i 16 

EN(KSSfIA)=EV(1 i IA)tVMU( 1 f IAf JKSSf ITF)+ 

.  EV(2»IA)tVNU(2>IAiJKSSiITF)+ 

.  EV(3.IA>*VHU<3fIAf  JKSSfITFH 

.  EV(4fIA)*VHU(4fIAfJKSSfITF)+ 

.  EV(5f  IAXVNUt 5>  I A»  JKSSr  ITF ) 

6  CONTINUE 

I  CONTINUE- 

DO  190  K=1 fNPNLU 
KSS=K+NPLH 
JKSS-LOC J(KSS ) 

DO  141  I A— 1 • 16 
KSSU=KSSUA( K > I A ) 

ENCKSSf IA)=XMA1  ( IAf JKSSf ITF) 4E(KSSf IA)+XHA2( IAf JKSSf ITF) tECKSSU 
1  fIA) 

141  CONTINUE 
190  CONTINUE 

C  2ND  HALF  OF  PROPAGATION  IN  LAX-UENDROFF  REGION 
DO  11  KSS=1 iNPLU 
JKSSzLOCJ(KSS) 

IF  ( ITF.LE.NLP<KSS) )  GO  TO  12 
DO  5  IA=2il6»2 
IIA=IREFA< IA) 

KSSU=KSSULH<K$SfIIA) 

E<KSSfIA)=PU1(IAfJKSSfITF)*EN(KSSfIA)+PU2(IAfJKSSfITF)»EN<KSSU 
1  >  IA) 

5  CONTINUE 
GO  TO  11 

12  DO  13  IA=2> 16» 2 
I IA=IREFA( IA) 

EV(1fIA)=EN(KSSfIA> 

DO  14  K=2f5 
LI=I+LXLI ( I IAfK) 

LJ= J+LXL J< I IAfK  ) 

KSL=L0CC<LIfLJ) 

EV(KfIA)=EN(KSLfIA> 

14  CONTINUE 

13  CONTINUE 

DO  16  IA=2r 16v2 

E(KSSfIA)=EV<1fIA)»PMU< If IAf JKSSf ITF) 

.  +EV<2fIA)*PNU(2fIAfJKSSfITF) 

.  +EV(3i IA) *PHU(3r IAf JKSSf ITF ) 

.  FEV( 4 1 1  A) tPHU( 4  f IAf JKSSf ITF ) 

.  +EV<5fIA)*PHU(5fIAfJKSSfITF) 

14  CONTINUE 

II  CONTINUE 

C  2ND  HALF  OF  PROPAGATION  IN  BOUNDARY  REGION 
DO  9  K=1fNPNLW 


KSS=K+NPLM 

JKSS=LOCJ(KSS) 

DO  41  IA=2. 16.2 
IIA=IREFA( IA) 

KSSU=KSSU A ( K » 1 1 A ) 

£(KSS»IA)=PU1(IA.JKSS»ITF)IEN(KSS . IA) +PU2( IA» JKSS. ITF )*EN(KSSU 
1  >IA) 

41  CONTINUE 

9  CONTINUE 

999  CONTINUE 
DO  7  KSS=l»NTOT 
DO  7  IA=2  r  16  f  2 
EN(KSS» IA)=E(KSS» IA) 

7  CONTINUE 

DO  6  I A  =  1 » 16 
DO  8  KSS=l.NTOT 

IF  (EN(KSS.IA).LT.O.)  EN(KSS. IA)=0. 

8  CONTINUE 
RETURN 
END 

FUNCTION  HBARF(E.UST) 

UST4=UST*»4 
EBAR=960400. BE/UST4 
HBARF=0.044*E  BARM  (-0.2) 

IF  (NBARF.LT. 0.008)  HBARF=0.008 

RETURN 

END 

SUBROUTINE  DOTPRD  (X.Y.Z) 

DIMENSION  X(16)>  Y ( 16  ) 

2=0. 

DO  1  1=1.16 
Z=Z+X( I >*Y< I ) 

1  CONTINUE 
RETURN 
END 


APPENDIX  C:  SAMPLE  WAVE  MODEL  RUN 


i 

1.  A  sample  flat-bed  wave  model  run  using  the  FORTRAN  listing  in 
Appendix  A  for  one  day  is  included.  The  input  consists  of  a  file  calle 
INPTEST  and  the  wind  file,  ILWND.  The  INPTEST  file  contains  all  the 
grid  information  and  the  parameters  and  constants  needed  to  run  the 
wave  model.  The  ILWND  file  contains  the  wind  speeds  and  directions 
needed  to  drive  the  wave  model. 

2.  The  output  included  in  this  report  contains  the  printout  of 
the  information  discussed  in  the  section  on  standard  output  from  the 
model.  Note  that  the  input  file  has  1-hr  time-steps  (TINC=3600.)  and 
has  wind  input  every  3  hr.  Stability  criteria  require  that  1-hr  time- 
steps  are  the  largest  that  can  be  used  with  this  input  setup.  The  out¬ 
put  prints  do  not  print  the  hourly  data--the  data  are  printed  every 

3  hr  to  match  the  input  wind  data. 

3.  File  INPTEST  describes  a  9x9  grid  with  a  grid  spacing  of 
222  km.  TINC  is  3600  sec  (1  hr).  The  frequency  range  is  valid  for 
oceanic  conditions.  One  special  output  station  at  1=5, J=8  is  included. 

4.  File  ILWND  is  the  input  wind  file.  The  72050100  is  the  date¬ 
time  which  means  May  1,  1972,  at  00  hours.  The  matrix  directly  below 
the  date-time  is  the  wind  speed  in  knots  and  below  this  is  the  direc¬ 
tion  matrix  (in  degrees)  that  corresponds  to  the  wind  speeds. 

5.  The  next  file  is  the  output  data  from  the  model.  The  prints 
have  been  labeled  to  make  this  information  self-explanatory;  but  if  any 
questions  should  arise,  see  the  section  on  standard  output  from  the 
model  on  page  16. 


File  INPTEST 


File  ILWND 


72050100 

18  21  20  16  10  8  10  8  8 

18  20  20  16  10  8  8  8  1 

20  21  20  18  14  8  8  10  1 

16  21  21  18  14  10  8  12  1 

12  21  20  fa  14  10  10  14  1 

14  21  18  16  14  12  10  12  1 

12  16  18  16  12  12  12  12  1 

6  14  16  18  14  14  14  12  1 

12  12  14  14  12  10  12  16  1 

155146137143142167203250250 
153143134141140165201248248 
145133126134149168200225225 
137125129128142164204232232 
122105122131139165195231231 
104107119125137178204232232 

50103123132140177197232232 

147119127135164188185218218 

332100124130154180188207207 

72050103 

15  21  22  18  10  8  10  9  9 

15  20  21  17  11  8  8  8  1 

16  19  20  17  14  8  8  10  1 

13  18  20  18  13  10  8  12  1 

11  17  19  18  13  9  9  14  1 

11  15  18  15  13  11  10  13  1 

11  12  17  14  11  11  11  12  1 

7  12  15  15  13  12  13  11  1 

10  11  13  13  11  10  11  15  1 

123139134143142169205254254 
121136131141140167203252252 
76122122131147170200229229 
212119125122141165202237237 
37  97118125141165193234234 
20105108120138177204234234 
337104111129141175201235235 
197112121137165188190224224 
316103131136158182194212212 
72050106 

12  21  25  20  10  8  10  10  10 

12  20  23  18  12  8  8  8  1 

17  16  20  18  14  8  8  10  1 

10  14  20  18  12  10  8  12  1 

10  12  18  18  12  8  8  14  1 

8  8  18  14  12  10  10  14  1 

10  8  16  12  10  10  10  12  1 

8  10  14  12  12  10  12  10  1 

8  10  12  12  10  10  10  14  1 

92132131142142171205257257 
90129128140140169203255255 

7110118128146171198231231 
285113122116141166199240240 
312  90114119142164190234234 
296104  98115139175203233233 
263105100126141174203235235 
245105115138165186194227227 
299106139141161182198215215 
72050109 


14  17  21  19  11  9  11  11  11 

13  15  20  17  12  9  9  9  1 

13  13  16  17  13  9  8  10  1 

10  10  15  16  11  9  8  10  1 

10  9  14  16  11  8  8  12  1 

10  8  13  13  11  10  10  12  1 

12  7  12  12  10  9  9  10  1 

10  8  11  12  12  10  10  8  1 

9  8  10  12  11  10  10  12  1 

65115125133135178210249249 
63113122131132176208247247 
1  73111121146181204233233 
313  491201 131 4S1 73202240240 
307179118124151172191232232 
287179106122143183201226226 
256173119132143179200225225 
246171138141166187192215215 
274167155144161179194206206 


72050112 

16 

12 

18 

18 

12 

10 

12 

12 

12 

14 

10 

16 

16 

12 

10 

10 

10 

1 

14 

10 

12 

16 

12 

10 

8 

10 

1 

10 

6 

10 

14 

10 

8 

8 

8 

1 

10 

6 

10 

14 

10 

8 

8 

10 

1 

12 

8 

8 

12 

10 

10 

10 

10 

1 

14 

6 

8 

12 

10 

8 

8 

8 

1 

12 

6 

8 

12 

12 

10 

8 

6 

1 

10 

6 

8 

12 

12 

10 

10 

10 

1 

38 

99119123127185212240240 

36 

961 

116121125183210238238 

355 

36103115146189207234234 

338345118110150180204238238 

301269123130160179191228228 

277254115129147189197216216 

247241139138145185196212212 

246238161143166186188201201 

247228171147161176189195195 

72050115 


19 

13 

14 

15 

9 

7 

10 

10 

10 

15 

11 

13 

14 

9 

7 

8 

9 

1 

14 

9 

8 

14 

11 

8 

7 

10 

1 

10 

7 

6 

11 

10 

8 

8. 

9 

1 

7 

8 

7 

11 

9 

7 

9 

12 

1 

8 

8 

6 

10 

10 

9 

10 

12 

1 

11 

7 

7 

10 

9 

8 

9 

10 

1 

12 

10 

9 

11 

12 

10 

9 

8 

1 

10 

9 

9 

11 

12 

10 

11 

11 

1 

46  75112134132182216254254 
44  73110132129179214252252 

15  29113126150185214231231 
7329156120147177212234234 

328270137136149183203229229 

271268178136144194208223223 

229245186145146189209220220 

228228190151167191204215215 

231225188158160183204207207 

72050118 

21  14  10  12  6  4  8  8  8 

16  12  10  12  6  4  6  8  1 


14  8  4  12  10  6  6  10  1 

10  8  2  8  10  8  8  10  ; 

4  10  4  8  8  6  10  14  1 

4  8  4  8  10  8  10  14  1 

8  8  6  8  8  8  10  12  1 

12  14  10  10  12  10  10  10  1 

10  12  10  10  12  10  12  12  1 

54  52106144136178218267267 
52  49103142134176216265265 
35  22122138155179218227227 
35311195131145174219228228 
354269151143138184213228228 
263280242143140197218227227 
210248231151146192221225225 
208216218158168194219227227 
213221203168159188217217217 
72050121 

24  16  8  10  7  5  6  6  6 

18  13  8  10  7  5  5  6  1 

14  74996581 

10  6  2  7  10  8  7  10  1 

4747968  14  1 

4  6  5  8  10  8  9  13  1 

7  7  7  9  9  9  10  12  1 

9  12  10  12  13  11  11  10  1 

8  12  9  11  13  11  12  12  1 

59  53  B71 27125153204274274 
57  50  84125122151202272272 
41  39  92130146154208224224 
42320192138139161211228228 
71255172148130166206230230 

226257226144134184217234234 

197236222150142181215233233 

196217216154165187211230230 

209224204163156182209217217 

72050200 

27  IB  6  8  8  6  4  4  4 

20  14  6  8  8  6  4  4  1 

14  64686461 

10  4  2  6  10  8  6  10  1 

4  4  4  6  10  6  6  14  1 

4  4  6  8  10  8  8  12  1 

6  6  8  10  10  10  10  12  1 

6  10  10  14  14  12  12  10  1 

6  12  8  12  14  12  12  12  1 

64  54  68109113128188280280 
62  51  65107111126186278278 
47  55  61122137129195219219 
49327188146133147202226226 
148239194154121148198230230 
188232208146128169214238238 
183222211149137168208238238 
183216213150162179201230230 
203226203157152174200214214 


Wave  Model  Output 


DATE-TIME=  72050100 

TIME  IN  HOURS  SINCE  START  OF  RUN  =  3.00 
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MEAN  DIRECTION  OF  THE  SEA(BEGREES)  =  225 
MEAN  DIRECTION  OF  ALL  WAVES  <DE0REES>=  231 
WIND  SPEED  AT  GRID  POINT(KNOTS)  =  14 
DIRECTION  OF  WIND  AT  GRID  POINT (DEGREES)1  231 


IiATE-TIflE=  72050100 
AL  GRlt'-SIGNlF ICANT  WAVE  HEIGHT  (HETEftE) 


o  o  o  o  o  o 


OOOOOOOO 


• 

• 

. 

• 

• 

• 

• 

• 

«-* 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

•  • 

o 

o 

o 

o 

o 

o 

o 

o 

o 

>-x 

O 

o 

o 

o 

o 

o  o 

o 

o 

o 

o 

o 

o 

o 

o 

O 

o  o 

<n 

UJ 

UJ 

o 

ro 

ro 

ro 

fl 

ro 

ro 

o 

ur. 

o 

ooinfi 

ci  rj 

00 

o 

o 

in 

in 

in 

in 

in 

in 

in  o 

* 

• 

• 

• 

• 

• 

• 

o 

• 

• 

• 

• 

• 

• 

• 

• 

• 

• 

o 

o 

o 

O 

o 

o 

o 

o 

o 

Ui 

o 

T 

CM  ro 

ro 

ro 

ro 

*-« 

o 

o 

-e 

♦ 

«r 

▼  o 

w 

cm  CM 

cj 

CJ  CM  CJ 

CJ 

ac 

o 

ro 

ro 

ro 

ro 

ro 

ro 

o 

o 

o 

in 

Cx 

m 

o 

o 

m 

m 

in  in 

in 

m  o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

►— 

o 

o 

o 

o 

©* 

o  o. 

00 

o 

o 

♦ 

♦ 

▼  o 

o 

CM 

CJ 

CJ 

*-« 

^  4 

Ui 

a- 

o 

ro 

ro 

ro 

ro 

ro 

ro 

▼ 

o 

cn 

m 

CD 

in 

CD 

Cx 

00 

o 

cn 

o 

UT 

UT 

in 

in 

m 

in  o 

• 

- 

• 

• 

• 

• 

• 

♦ 

• 

• 

• 

• 

• 

• 

• 

• 

eu 

• 

• 

• 

• 

•  . 

o 

o 

O 

o 

O 

O 

o 

o 

o 

U- 

o 

-o 

-o 

'O 

-o 

Cx 

Cx 

CD 

o 

z 

o 

♦ 

^  o 

o 

o 

(_> 

z 

UI 

o 

ro 

V 

ro 

o 

o 

o 

o 

Ox 

CJ 

o 

Cx 

o 

o 

CO 

o 

m 

in 

bT 

in 

u"> 

in 

in  o 

o 

O 

o 

o 

o 

o 

o 

o 

o 

►— 

o 

«r 

ro 

ro 

'O 

o 

1=) 

o 

▼ 

«r 

♦  o 

CJ 

•-« 

o 

UJ 

*— « 

l-t 

a: 

o 

♦ 

o 

M 

o 

«r 

CO 

m 

CJ 

UT 

o 

UJ 

o 

b-) 

in 

in 

in 

b”) 

in 

m  o 

i=a 

• 

• 

u. 

o 

o 

o 

o 

o 

o 

o 

o 

o 

ro 

CJ 

ro 

CJ 

ro 

ro 

o 

o 

•*- 

▼ 

o-  o 

z 

« 

<r 

UJ 

UJ 

o 

i n 

m 

m 

in 

♦ 

o 

z 

o 

>o 

o- 

CJ 

o- 

ro 

o 

Cl. 

o 

in 

in 

m 

in 

in 

m 

in  o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

ro 

CM 

CJ 

CJ 

•»j 

CJ 

CJ 

o 

A 

o 

♦ 

o 

•M 

»-4 

•j 

ur 

cc 

o 

a 

o 

in 

n 

in 

in 

bO 

o 

o 

ro 

ro 

bT 

in 

Cx 

ro 

o* 

o 

o 

O- 

O' 

cc 

O' 

Ox 

Ox 

Ox  o 

_ j 

• 

—I 

o 

o 

o 

o 

o 

o 

o 

o 

o 

<z 

o 

ro 

CJ 

o 

o 

O 

o 

<c 

o 

m 

m 

m 

in 

b-> 

bO 

in  o 

M 

K— 

■  4 

•“* 

«*j 

M 

1— 

<r 

<x 

o 

o 

o 

o 

o 

o 

o 

u. 

o 

o 

o 

o 

o 

o 

U- 

o 

o 

o 

o 

O 

o  o 

DATE-  f  IH£  =  72050103 

Fine  IN  HOURS  SINCE  START  OF  RUN  =  6.00 
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MEAN  DIRECTION  OF  ALL  WAVES  >DEGREES)=  23J 
WIND  SPEED  AT  GRID  POINTIKNOTS)  =  14 
DIRECTION  OF  WIND  AT  GRID  POINT ( DEGREES >=  23J 
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DATE-TIME*  72050106 

TINE  IN  HOURS  SINCE  START  OF  RUN  =  9.00 
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APPENDIX  D:  NOTATION 


Boundary  point  in  grid 
Group  velocity 

Grid  spacing  in  wave  model  (x-direction) 

Spacing  between  grid  points  in  the  longitudinal  direction 
(only  when  using  a  spherical  orthogonal  grid) 

Energy  associated  with  wave 

Average  nondimensional  energy  at  a  specific  point 
Frequency  of  wave 
Frequency  of  spectral  peak 

Minimum  frequency  that  is  input  into  wave  model 

2 

Gravitational  acceleration  (980  cm/sec  ) 

Significant  wave  height 

Average  wave  height  associated  with  E 
Latitude  counter 

Latitude  counter  of  center  of  grid 
Lax-Wendroff  water  point  in  grid 
Number  of  columns  in  grid 
Number  of  rows  in  grid 
Number  of  grid  points 
Grid  spacing  in  degrees 
One-dimensional  energy 
Time 

Nondimensional  time 
Time  increment  in  wave  model 
Wind  speed 
Friction  velocity 

Difference  between  two  separate  frequencies 
Direction  of  wave 
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1.  Page  A6 ,  after  line  38  (402  DO  20  KSS=1 ,NPTSP1) ,  add: 

SCAA  (KSS)=0 . 0 
SWLSC  (KSS)=0 . 0 

2.  Page  A10,  after  line  18  (DO  66  KSS=1,NPTS),  delete: 

HI GHE=HSCALE (KSS ) 

3.  Page  All,  after  line  39  (123  KFRQ=AF07(KSS) ,  add: 

H I GHE=HSCALE (KSS ) 

4.  Page  A15,  line  26  (IF(LA(I ,J) .EQ.O)  L0CC(I , J)=NPTSP1) ,  change: 
.EQ.  to  .LE. 

5.  Page  A15,  line  55,  change  VMU(3,IA,ITF)=EMU2*B*EMUP  to: 
VMU(3,IA, ITF)=-EMU2*B*EMUP 

6.  Page  A16,  line  11,  change  PMU(3 , ITF)=EMU2*B*EMUP  to: 

PMU(3, ITF)=-EMU2*B*EMUP 
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